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Abstract 


This  project  designed  nanoscale  solar  voltaic  cells  by  using  carbon  nanotubes.  The  Field- 
Effect  Transistor  of  canbon  nanotube  was  adopted  and  asymmetric  Schottky  barriers  made 
of  two  different  metal-semiconductor  barriers  at  both  ends  of  carbon  nanotubes  were 
employed  to  enhance  the  Einstein  photoelectric  electron-hole  separation.  Two  architectures 
of  the  nanoscale  solar  cell  were  designed  and  analyzed.  This  project  also  examined  the  use 
of  finite  element  analysis  to  simulate  the  process  of  nanoindentation  to  gain  insights  into 
the  behavior  of  elastoplastic  and  viscoelastic  materials.  The  study  examined  the  behavior  of 
a  calibrated  elastoplastic  material  model  during  indentations  using  different  indenter  tip 
geometries.  The  study  then  compared  simulations  and  experimental  results  for  the 
elastoplastic  material.  Additionally,  the  study  examined  the  calibration  of  a  viscoelastic 
material  model  from  reported  data  and  its  ability  to  predict  indentation  response. 
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I.  Summary  of  work  in  nano  scale  solar  voltaic  cells  by  using  carbon  nano  tubes 

1.1  Introduction 

Shockley  and  Queisser  have  computed  in  1961  the  upper  limit  of  30%  converting 
1000  Watts/nr  sun  to  300  Watts/m“  of  any  Schottky  semiconductor,  e.g.,  silicon  p-n 
junction  device  at  the  band  gap  at  1.1  eV  as  discovered  by  Bell  Lab  in  1954;  and  such  a 
solar  battery,  installed  by  Hoffman  electronics,  on  board  of  US  Satellite  Vanguard  1,  1958. 
The  upper  limit  is  due  to  photo-excitation  mechanism  which  must  be  above  the  Fenni  level 
of  the  silicon  or  the  GaAs  semiconductor  band  gap  at  1.1  eV  or  1.52  eV  respectively,  which 
excludes  the  long  infrared,  microwaves,  radio,  except  visible  photons  to  generate  free 
conduction  band  electrons.  In  this  project,  we  explored  Carbon  NanoTubes  (CNT) 
technology  to  break  the  Shockley-Queisser  theoretical  barriers  of  visible-only  solar  cells. 
We  adopted  the  simple  spherical  camera  lens  which  is  made  of  Pb-Crown  glass  of  F#=0.7 
af  the  Lensdiameter  =2cm,  which  will  provide  a  unifonn  focusing  of  sunlight,  from  the  visible 
to  near  infrared  (NIR),  up  to  millions  times  stronger  upon  at  the  yellow-light  focal  plane, 
enjoying  the  maximum-number  of  photons. 

1.2  Design  of  voltaic  solar  cells  using  carbon  nanotubes 

The  diameter  of  lenslet  array  was  chosen  as  2  cm  and  the  lenslet  was  made  of  Pb-lead- 
Crown  glass  with  the  refractive  index  around  1.52  and  F#  is  0.7,  focusing  the  maximum 
number  of  photons  in  the  visible  yellow  0.635  pm  and  red  0.73  pm,  at  the  sharpest 
Rayleigh  diffraction- limit  spot  size  yellow-light  focal  plane.  Thus,  judging  from  the  Planck 
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radiation  spectrum  shape  equivalent  to  the  sun  at  6000  °K  (peak  at  Yellow)  -5000  °K  (peak 
at  Red),  without  the  atmosphere  in  the  space,  the  visible  band  was  considered  from  0.4pm 
to  0.73pm,  and  the  infrared  band  from  0.73  pm  to  1.12  pm  being  limited  by  the  Pb-Crown 
glass  flat  transmission.  Then,  thermal  re-remittance  loss  due  to  a  hot  solar  panel  had  been 
considered  to  be  negligibly  small  percentage  by  Landsberg  and  Baruch.  Furthermore,  we 
realized  the  band  gap  concept  of  semiconductor  to  be  a  lower  threshold  of  minimum 
photon  energy  to  excite  conduction  charge  carriers,  and  thus  we  chose  a  relatively  larger 
cross-section  and  lower  band-gap  threshold  active  element.  It  turns  out  to  be  the  near 
infrared  at  1.12  pm  and  Carbon  NanoTube  (CNT)  with  the  diameter  of  0.66  nm.  The  focal 
plane  was  chosen  for  Yellow-light  at  the  Rayleigh  diffraction  limit  spot  size,  which  had  the 
maximum  number  of  photons  concentrated  at  the  focal  plane. 

Mosaic  solar  panel  at  the  size  of  2m  by  2m  could  be  made  manageable  to  a  Space- 
Robotic,  e.g.  Canadian  ann.  It  was  packed  with  100x100  unit  volume  cells,  as  shown  in 
Figure  1. 


F#  =  Z/Diens  =  0.7 

Aiens  =  7rr2=  n  (0.01m)2-=  0.0003 14m2=  n  cm2 

honev  -  comb  lenselet  collector 

active  bill  b  actor Panel  = - : - 

Squared  Frame  Panel 

_  NxMx/r-r2  _  NxMxtt 

[( r(N  -  1)V3  +  2 r)  •  (2 M  + 1  )r]  [(^3 N  +  2-S)-  (2 M  + 1)] 


(1) 

(2) 


(3) 


4 


When  N  and  M  become  very  large,  the  solar- Active  Factor  panel  is  near  a  constant  of 
7i/2*3A  «  90%.  Thus,  the  total  solar  energy  is  at  most  reduced  by  10%  frame  boundary 
used  to  hold  the  lenslet  array  in  to  the  firm  place. 


Lenslet 


Panel 


V  Mxcells 

=  (2A/  +  I>- 


=  r(N-\)S+2r 


Figure  1 .  Solar  cells  fill  in  a  percentage  of  the  panel  Crown  glass  is  type  of  optical 
glass  used  in  lenses  and  other  optical  components.  Crown  glass  is  produced  from 
alkali-lime  (RCH)  silicates  containing  refractive  index  (-1.52)  and  low  dispersion 


In  order  to  estimate  the  relative  output  electric  power  based  on  our  several  design 
architectures,  we  could  fix  the  common  microscopic  parameters,  such  as  the  photoelectric 
surface  work  potential  function  WEimtem  and  the  difference  effect  of  asymmetric  potential 
barriers  d>schottky  at  both  ends,  as  well  as  the  Einstein  photoelectric  conversion  efficiency, 
p Einstein  of  CNTs  and  the  carrier  transport  collision  extinction  ExcomSion(to  be  optimized  by 
our  subsequent  kinetic  R/D).  In  a  unit  cell  system,  the  focus  plane  was  generally  chosen  at 
the  visible  yellow  light  in  favor  of  generating  a  higher  velocity  Vvis  of  the  electrical  charge 
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carrier  which  contribute  to  a  larger  photo-electric  current  using  the  justifiable  classical  and 
semi-classical  formulae: 


Iv»  =q  Ne  Vvis.  (4) 

'ta‘  “  ~  ^  ‘I^Sbookley  “  ^^vts  (5) 

fellow  =  8.2x  105m/s,  Vnir  =  6.2x  105m/s  (6) 

In  our  design  we  approximated  the  integrated  majority  yellow  photon  visible  band  at  a 
typical  yellow  light  and  the  Crown-glass-passable  near  infrared  light  NIR  as  follows: 

^yellow  —0.635pm,  Eyeiiow  =  1.26xl02  J/(  s-m2),  Oyeiiow^  2.97x  1015rad/s  (7) 

and  for  NIR  light 

Xir  =1.12pm,  Eir  =  78  J/(s  m2),  coffl  =  1.68x  1015  rad/s  (8) 

A  more  photon  collector  in  a  larger  CNT  diameter  associated  with  infrared  (IR)  photons 
according  to  Seito  band  gap  relation: 

E gaitQ  Bq  —  —  —  ,  EEinstein  hcA,=(1.24  eV)/(X  pm)  (9) 

a  CNT 
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where  ac-c  is  1.41  A,  neighborhood  hoping  distance,  yo  =2.6  eV=  energy  of  neighborhood 
hopping  interaction. 

We  chose 

dCNT(^NiR)=  0.66  nm;  EBG(NIR)=1.107  eV  (10) 


The  band  gap  formula  in  Equation  9  was  verifiable  because  of  the  transitional  regime 
permitted  a  mixed  treatment  of  classical  physics  and  quantum  mechanics.  The  propagating 
electron  along  CNT  whose  quantum  mechanical  de  Broglie  wavelength  going  around  the 
circumference  of  CNT  was  computed  as  Sommerfeld-Bohr  quantization  condition  as  an 
integer  of  the  circumference  as  follows: 


**"■*"'  -  rJmm  - 

n  7f  dc.Yr(A.Ytf?) 


n=  1,2,  etc.  (11) 


de  Broglie  matter  wave  going  around  single  layer  of  CNT  lattice  equals  once  to  be  in  phase 
with  itself  for  n=l  for  NIR,  and  then  for  n=2  twice  for  the  visible  red. 

Because  of  a  larger  basket  can  catch  both  the  smaller  visible  photons  generating  larger 
kinetic  velocity  of  the  charge  carriers  and  the  larger  NIR  photons  generating  smaller  kinetic 
velocity  of  the  charge  carriers,  according  to  the  semiconductor  Femi  level  at  NIR  band  gap. 
However,  whenever  possible,  we  chose  the  focal  length  to  be  located  at  the  Lord  Rayleigh 
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diffraction  limit  smallest  spot  size  at  the  shorter  yellow  wavelength  where  an  abundance  of 
solar  photons  existed,  and  produces  larger  velocity  charge  carriers  than  NIR. 

We  suppose  that  all  16  NIR-CNTs  may  be  spirally  packaged  in  the  Rayleigh  yellow-light 
focusing  origin,  of  which  the  central  core  thickness  becomes  increased  linearly  from 
0.66nm  to  0. 66x16=10. 56nm,  where  the  diffraction-limited  spot  size  was  de-focused  from 
for  NIR. 


Given  the  common  material  constants  in  MKS  units  adopted  by  solar  energy  community: 


me  =  9.1  lxlO'31  Kg,  h  =  1.054xlCr34  J-s,^  =  -1.602  x  10'19Coulomb  (12) 


We  computed  the  total  number  of  photons  involved  in  photoelectric  conversion  by  the 
following  conservation  of  energy  formula 


Total  _  Energy  {shorter  &  including  A) 
=  Ei-  Au,v  ■  -p22-  =  Nplmlm  ha 

SpotSize 


~\T  _  T y  a  "^CNTs  (  ^yellow )  jo- 

^Vphoton  yellow-  & yeuow  *  ALens  '  77  n  ’  W yellow 

A  SpotSize  V  ^ yellow  ) 

=  5.44  xlO16 


JVpto,„„MR=  /  /j.raM~1.89xl01! 

^SpotSize  \  ^NIr) 


(13) 


(14) 

(15) 
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We  have  the  larger  IR-diffraction  spot  size  and  larger  photon  collector 


Dspotsize(^NiR)  -0.45pm,  dcNi(^NiR)-  0.66nm,  (16) 

as  opposed  to  the  smaller  yellow  diffraction  spot  size  and  smaller  photon  collector  at  the 
visible  yellow  band  gap  CNT : 

DspotSizc(/-ycii o w)  —  0.051pm,  d(\i( /'.yellow)  0 . 3 7 5 n m  (17) 

1.3  Two  architectures 

We  compared  the  relative  conversion  efficiency  of  one-layer  architecture  with  that  multiple 
layers,  bearing  in  mind  but  postponed  until  the  final  tradeoff  consideration  between 
manufactory  difficulty  and  the  cost. 

(i)One  Basic  Layer  Architecture: 

The  architecture  contained  16  NIR  CNT  at  0.66  mn  diameter  is  shown  in  Figure  2. 


One  layer  system 


Figure  2.  Based  on  the  architecture  picture,  we  computed  the  total  cross  section  area  that 
the  visible  yellow  light  was  intercepted  by  the  NIR  CNT  in  a  spiral  stack  of  16  on  the 
visible  yellow  focal  plane. 
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d cntj  NP  ■  dCNT  (A)  •  DSpotSbe  aoverlap  j  ■  (N p  1) 


=  A 


BG 


0-9/n  •  ac 


■NP~ 


(2r0-^c)2 

1.54 


(NP- 1) 


^cnts(^nir-1-12  ju.ni)  —  0.00476pm^, 

^CNTs  (^yeiiow=0.635|jm)=0.0027jjm2  (18) 


where  the  Rayleigh  spot  size  was  focused  at  yellow  light  focal  plane,  which  was  smaller 
than  that  of  NIR  as  indicated  on  the  first  layer. 


^SpotSize  (^-yellow)  0.051pm*",  A spotSizeO'-NiR ) — 0.318pm 


(19) 


We  considered  a  detail  kinetic  theory  of  the  charge  migration  suffering  the  collision 
extinction  with  the  neighborhood  carrier  pairs. 


N 


phone 


■Tim 


'  (1  EX  CoUjsion  )  -  Nt 


I  Total  TfN carrier  NIR  '  V^IR  qN carrier  yellow  ’ 

/totalonelayer  at  Yellow  FP/  ^1  Einstein' {,  1 -Excoujsion)  9.22x10  Amp.  Ill  MKS  (20) 


10 


One  tayef  r/ttam 


'bp’rim 


MMytf 


Thickness 

=8dcNT 


Figure  3.  The  fan-out  pattern  of  our  volume  cell  is  horizontally  along  the  radial  direction 
at  every  22°  each,  the  closed  packaging  of  Np=16  NIR  CNTs  must  be  overlaid  passing 
one  another  in  the  central  core  creating  a  stack  of  tubes  like  a  fireman  spiral  staircase. 


(ii)  Two  Layers  Architecture 

We  assumed  the  kinetic  velocity  of  a  single  charge  carrier  to  be  identical  at  either  layers 
and  its  velocity  depends  sole  on  a  single  photon  energy  which  has  excited  the  carrier  out 
of  the  Fermi  level  after  overcome  the  Einstein  surface  work  potential  and  transporting 
along  the  net  field  of  two  asymmetric  Schottky  barriers  at  both  ends:  Palladium  (Pd)  as  P- 
type  and  Aluminum  A1  as  N-type. 


We  focused  the  yellow  light  at  the  focal  plane  where  8  NIR  CNTs  were  laid  in  the  first 
layer,  then,  the  rest  near  infrared  (MIR)  light  was  focused  at  the  second  layer  focal  plane 
where  8  NIR  CNTs  were  laid. 
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2  meUVIS  ~  fo^VIS  qW Einstein  +  ^  Schottky  —  VIS 


Layer  1 


Electrode 


Layer  2 


Visible  IR 


2  CNT  layers 


Figure  4.  Design  model  of  3D  lens  cavity  consisting  of  8x2=16  CNTs. 


As  shown  in  Figure  5,  the  visible  light  was  focused  on  the  first  yellow  focal  plane  laid  with 
8  NIR  CNTs  plane 


A cnjts  -  0.00134pm2,  4spotSizc~  0.05 1  pur 


(21) 


TWo  layer  system 


Figure  5  Two  Layer  Architecture  indicated  4  focal  spot  sizes  covered  by  8  NIR  CNT 
photon  collectors 
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NIR  light  was  defocused  on  the  first  yellow  focal  plane  laid  with  8  NIR  CNTs  plane 


yfcNT  =  0.00334pm2,  ^sPotsize=0.3 16pm2 


(22) 


Yellow  light  became  defocused  on  the  second  NIR  focal  plane  laid  with  8  NIR  CNTs  plane 


4  (NT  -0.  00 1 9pm2,  zlSpOtsize=0. 1 02pm2 


(23) 


NIR  light  was  focused  on  the  second  NIR  focal  plane  laid  with  8  NIR  CNTs  plane 


4cnt  -0.0024pm2,  z4spotsize— 0. 1 58pirr 


(24) 


Then  the  total  energy  had  two  band  widths,  yellow  and  NIR,  spreading  at  the  first  yellow- 
light-focal  plane,  produced  different  numbers  of  photons: 


Total  _  Energy  llow  _  at  V  _  Layer 


=  E 


yellow  _  1  st  _  layer 


(4.M 


A 


CNT  yellow 


Lens 


A 


=  N 


photon  _  yellow 


I  -ti-  CO 


(25) 


yellow 


SpotSize  _  yellow 


Total _ Energy N1R  _at  f'  _layer 


=  E 


NIR  _  1  st  _  layer 


(A NIR  )  ‘  A 


A, 


CNT  1R 


Lens 


A 


~  N photon  _IR_  1  ’  h  '  <®nir 


(26) 


SpotSize  _  IR 
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The  remaining  non-intercepted  spectrum  energy  is  re-distributed  in  the  second  NIR-focal 
plane 


Total  _  Energy yellow  _  2nd  _  layer 


=  E 


yellow  _  2nd  _  layer  \  yellow 


^  A  CNT  _  yellow  _  AJ- 

\A yellow)'  A  Lens  ,  ~  ^  > 

, SpotSize  _  yellow 


photn  _  yellow  _  2  ^  ^yellow 


(27) 


Total _Energy NIR  _2"‘  _layer 
=  E,b  ' . -  (2MTR )  •  ALens  ■  — 


IR  2nd  _  layer 


CNT  NIR 


N ph0ton  NIR _2  '  ^  ‘  ^NIR 


SpotSize  _  NIR 


(28) 


The  total  number  of  charge  carriers  may  be  estimated  in  terms  of  the  common  efficiency 
parameter  and  the  transport  collision  extinction  probability. 


N. 


carrier  _  yellow 


=  (N 


+  N 

photon  _  yellow  _  1  photon  _  yellow  _  2 


)  collision  )  Einstein 


(29) 


N, 


carrier  NIR 


=  (7V 


photon  _  NIR 


i+tf, 


photon  _  NIR  _  2 


)  •  (1  Excollision  )t)Ein 


(30) 


I  =a-(N  -V  +N  ■ V  ) 

total  two  layer  j  V  carrier  NIR  NIR  carrier  yellow  yellow  I 

N photon_yellow_l  2.7x10  ,  N photonNIR  1=  1.39xl016, 

=>  fVphoton_yellow_2==  1 . 8 6 X  1 0 * Aphoton  NIR _2=  1  -96x  1  O'6  (31) 


ftotal  two  layer/  Einstein  ( 1  -Excoiiision)  =  9.42x  1 03  AMP  in  MKS  (32) 


14 


1.4  Discussions  and  conclusions 


The  physics  of  Diffraction  limit  of  focus  spot  size  is  estimated  by  Lord  Rayleigh  about 
the  diameter  length  R(z)  ,  Z0  is  the  focus  length 


when  Z=Zn 


R(z)  = 


fZ-Z03 

1 

2 

( z ) 

2 

1  + 

— 

l  n  ) 

Uo  J 

R0  = 


2A-Z0y 
V  n  J 


(33) 


(34) 


It  could  be  expressed  to: 


0.5  kX  X 

NA  /#<2  ~  2 


(35) 


where  k  is  a  constant  of  spot  size,  NA  is  numeric  aperture  (n -sin  6).  F#  is  of  F  number 
=(Zo/D),  D  is  lens  diameter.  We  can  reduce  the  spot  size  R  by  tailoring  the  refractive 
index  n  and  increasing  the  lens  diameter,  implying  a  possibility  of  magnification  of  the 
visible  solar  light  within  a  smaller  spot  size  than  that  of  IR. 


15 


Based  on  the  Lord  Rayleigh’s  diffraction-limited  focal  length  known  as  the  Rayleigh  range 


x  (Spot  size  radius  defocus  becoming 


V2 


factor  distance  indicated  in  Figure.  6). 


Figure  6.  Lord  Rayleigh  Diffraction  Limit  Focal  Spot 


Thus,  we  can  find  the  wavelength  at  0.73  pin,  the  Rayleigh  range  x=361  mn,  when  the 
wavelength  is  1.12|um,  the  Rayleigh  range  x=  568  mn.  Therefore,  the  thickness  of  a  stack 
of  spiral  NIR  CNTs  has  the  total  thickness:  16  x  0.66nm=l  0.56  mn  which  is  about  1%  of 
Rayleigh  range.  Consequently,  the  spot  size  variation  of  each  NIR  CNT  along  the  same 
layer  may  be  neglected. 


16 


One  layer  system  puts  16  larger  photons  collectors  made  of  NIR  CNTs  onto  the  yellow 
light  focal  plane,  compared  better  than  that  of  two  layers  having  8  NIR  CNTs  on  each 
layers. 

Antal  two  layer  / 7total_one_layer“  =  1.02  (37) 

To  increase  the  efficiency,  we  exploited  the  transition  regime  matter  of  which  the  Carbon 
NanoTube  (CNT)  is  a  typical  example:  the  charge  carriers  along  the  elongated  axial 
direction  follow  the  classical  physics  model,  while  the  azimuth  cross  section  of  CNT  is 
governed  by  the  quantum  physics  lattice  vibrations,  phonons.  A  higher  energy  photon 
which  can  pass  the  CNT  threshold  of  higher  band  gap  and  producing  higher  velocity 
electric  current  j,  along  the  longitudinal  direction  which  as  we  said  before  is  valid  in 
classical  kinetic  theory  of  current,  j=qpv,  which  is  proportional  to  the  charge  q,  the 
density  p,  and  the  velocity  v.  However,  the  quantum  physics  is  also  valid  around  the  CNT 
tube  cross  section,  where  the  de  Broglie  wavelength  of  electron  wave  function  is  defined 
as  p =h/X.  We  must  demand  the  uniqueness  of  the  wave-function  when  it  is  going  around 
the  circumference  of  a  CNT  tube,  which  is 


/l 


'DeBoglie 


electron 


n ' 7i ' d cnt (Amr ),n  1 ,2,  etc. 


(38) 


as  Sommerfeld  quantization  rule  has  done  for  Bohr  Hydrogen  atom,  and  we  have  chosen 
n=l  for  the  ground  state  introducing  the  diameter  d  of  CNT.  Thus,  we  have  verified 
independently  Saito  band  gap  formula 

EsG~l/dcNT-  (39) 
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This  argument  demonstrated  the  nanotechnology  of  CNT  to  be  valid  in  the  transition 
regime  of  both  classical  and  quantum  mechanical  physics.  Thus,  we  can  design  from  the 
microscopic  semi-conductor  Einstein  photoelectrical  effect  to  the  macroscopically 
demonstrated  Photovoltaic  Solar  Cell  (PVSC)  demonstrated  in  Carbon  NanoTube  (CNT) 
first  by  Li  of  GE  with  2-5%  efficiency.  We  took  advantage  four  more  considerations 
beyond  GE  experiment:  the  full  spectral,  the  full  polarization,  the  focusing  lens,  andthe 
asymmetric  Schottky  effect. 

The  solar  cells  enjoy  seven  attributes:  (1)  the  design  of  multiple  spectral  layered 
architecture  under  the  unit  optical  axis  of  one  simple  focusing  lens  was  based  on  (i)  GE, 
Li  APL  2004,  demonstrated  CNT  as  PVSE  with  5  V  dc-bias  and  2  symmetric  Schottky 
barriers  read-outs,  2%~5%  quantum  efficiency  (QE)  for  mono  or  full  polarization;  (ii) 
Saito  computed  the  CNT  band  gap  Eg- 1 /Diameter.  (2)  We  adopted  the  Field-Effect 
Transistor  of  CNT,  without  the  undesirable  dc-bias  supply  battery  in  the  Space;  instead, 
we  employed  asymmetric  Schottky  barriers  made  of  two  different  metal-semiconductor 
barriers  at  both  ends  of  CNT  to  enhance  the  Einstein  photo-electric  e-  •  -h+  separation  to 

achieve  the  CCD-like  dual  read-out.  (3)  We  exploited  the  Rayleigh  diffraction  limit  focus 
spot  size  as  a  function  of  the  spectral  wavelength  in  Np  cross-polarization-bars.  (4)  over 
full  spectral  bands  in  Ns  layers  avoiding  the  shadowing  effect  to  enhance  the  total  QE 
without  the  wasteful  fill  factor.  (5)  1-D  CNT  enjoys  less  thermal  noise  in  wasted  heat 
-1/2  KBT=T/40  eV.  (6)  ~0(N)  PIN-OUT  dual  CCD-like  Circuitry.  (7)  We  minimize  the 


kinetic  recombination  of  neighborhood  electrons  and  holes  in  the  phase  space  of  the  full 
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2D  surface  with  spiraling  charge  carriers.  Our  analysis  of  the  optics  throughput  efficiency 
indicated  the  trade-off  design  of  a  circuit  of  2  different  metal-semi-conductor  Schottky 
Barrier  contact  electrodes  for  16  CNT  =  Ns(2)  •  Np(8)  sources,  rather  than  original 

proposal  of  16  CNT=4x4  in  the  3D  cavity  structure  which  tracks  the  bare  sun 
automatically  in  the  Space  and  absorbs  the  dominated  solar  spectrum  energy  and  thus  can 
increase  the  total  quantum  efficiency. 
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II.  Summary  of  work  in  nanoindentation  of  elstoplastic  and  viscoelastic  materials 

II.  1  Nanoindentation  in  an  Elastoplastic  Material 

Nanoindentation  is  an  experimental  method  used  to  characterize  materials  at  the  nano 
and  micro  scale.  During  a  nanoindentation  test  a  diamond  indenter  tip  is  pressed  into  a 
material  up  to  a  specified  force  or  depth  and  then  withdrawn.  As  the  test  is  performed,  the 
indentation  force  and  depth  are  recorded  and  results  are  normally  presented  in  the  form  of  a 
load-displacement  curve.  This  load-displacement  curve  may  be  used  to  detennine  different 
material  properties  and  characteristics  (Oliver,  2004;  Mencik,  2007).  An  important  factor  in 
examining  the  results  is  the  type  of  tip  used  in  indentation.  Berkovich  indenters  are  one  of 
the  most  common  types  of  indenter  tips  used.  A  Berkovich  tip  is  faceted  in  the  form  of  a 
three  sided  pyramid.  Information  relating  to  stress-strain  behavior  in  an  elastoplastic 
material  during  loading  and  unloading  may  be  provided  by  finite  element  analysis.  The 
commonly  used  Berkovich  indenter  with  its  three  sided  pyramidal  geometry  must  be 
modeled  and  analyzed  as  a  3D  contact  problem  in  finite  element  simulations.  However,  the 
common  approach  in  finite  element  simulation  is  to  use  a  2D  axisymmetric  indentation 
model  which  uses  a  conical  shaped  indenter  as  an  equivalent  to  a  Berkovich  indenter 
(Knapp,  1999;  Panich,  2004;  Pelletier,  2000).  The  main  reason  for  this  approximation  is 
that  the  3D  contact  problem  may  require  a  large  amount  of  computational  resources  and 
usually  takes  a  long  time  to  complete,  especially  when  using  a  fine  3D  mesh  for  accurate 
results.  To  overcome  these  challenges,  a  70.3°  half-angle  conical  indenter  is  often  used  as  a 
substitute  for  modeling  the  Berkovich  tip,  on  the  bases  that  it  provides  the  same  indentation 
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depth  to  projected  area  ratio  of  a  Berkovich  indenter.  By  using  a  conical  indenter  the  finite 
element  analysis  can  be  simplified  by  using  a  2D  axisymmetric  model  which  allows  for 
quicker  simulations  using  fewer  computer  resources.  Since  this  is  a  common  approach  for 
the  finite  element  analysis  of  indentation  problems  it  is  important  to  compare  the  finite 
element  analysis  results  for  both  conical  and  Berkovich  indentations  to  investigate  any 
deviations  in  behavior  between  the  two.  While  there  have  been  several  recent  studies 
comparing  finite  element  simulations  for  the  two  indenters,  the  results  presented  by 
different  investigators  are  not  in  full  agreement. 

The  objective  of  this  thesis  is  to  present  a  detailed  finite  element  study  of 
nanoindentation  in  order  to  gain  additional  insights  into  the  effect  of  tip  geometry  on  the 
force-displacement  curves  obtained  in  nanoindentation  tests.  The  previous  studies  used 
different  elastoplastic  models  including  the  use  of  different  power-law  hardening  and  the 
use  of  elastic-perfect  plasticity.  This  study  will  investigate  the  behavior  of  conical  and 
Berkovich  indenters  for  materials  obeying  a  combined  isotropic  and  kinematic  hardening. 

II.  1 .2  Literature  Review 

One  of  the  earliest  studies  that  numerically  compared  different  indenter  tips  was 
perfonned  by  Li  et  al.  (2004).  In  their  study  they  compared  several  types  of  indenters, 
Berkovich,  Vickers,  and  Knoop  to  several  different  conical  indenters  with  half  angles  that 
were  determined  by  different  equivalency  rules.  The  study  provides  simulations  with 
indentation  depth  of  approximately  lOOOnm  in  an  elastoplastic  material  with  isotropic 
hardening.  The  load-displacement  curves  by  Li  et  al.  (2004)  show  close  agreement  between 
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the  Berkovich  result  and  the  conical  result  with  the  Berkovich  load  being  slightly  higher  in 
the  mid  range  of  indentation  depth.  The  authors  conclude  that  stress  fields  between 
Berkovich  and  conical  indenters  are  different  and  if  stress  field  is  of  interest  a  3D 
simulation  is  called  for.  However,  if  only  load-displacement  curve  is  needed  then  the 
Berkovich  indenter  can  be  sufficiently  simulated  by  a  conical  indenter  (Li,  2004). 

Swaddiwudhipong  et  al.  (2006)  also  attempt  to  address  equivalency  of  Berkovich  and 
conical  indenters  by  performing  2D  and  3D  finite  element  analysis  to  simulate  the  load- 
displacement  response  in  elastoplastic  materials.  Indentation  simulations  were  performed 
for  both  conical  and  Berkovich  indenters  using  different  materials,  including  6061 
aluminum,  that  are  set  to  obey  power-law  strain  hardening.  Indentations  were  simulated  to 
a  relatively  deep  5pm.  The  paper  provides  a  load-displacement  curve  comparing  conical 
and  Berkovich  indentations  in  6061  aluminum.  In  this  curve  the  Berkovich  indentation 
curve  is  noticeably  lower  than  that  of  the  conical  indentation.  The  authors  conclude  that  the 
equivalency  between  the  two  types  of  indenters  in  terms  of  the  curvature  of  the  load- 
displacement  curve  is  not  valid  across  the  range  of  material  properties  considered 
(Swaddiwudhipong,  2006). 

Finite  element  simulations  of  nano  indentation  in  elastoplastic  materials  with  no 
hardening  were  performed  by  Xu  and  Li  (2008)  to  study  the  effects  of  indenter  geometry 
and  mechanical  properties  of  the  indented  material.  The  indenter  geometries  considered 
were  conical  and  Berkovich.  The  study  was  specifically  looking  at  correction  factor  but 
there  are  also  plots  that  directly  compare  load-displacement  behavior  of  both  indenter 
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types.  In  the  provided  plots  the  Berkovich  indentation  curves  are  noticeably  higher  than 
those  for  the  conical  indenter  across  the  range  of  material  properties  shown  in  the  paper. 
The  difference  is  explained  as  a  geometry  effect  with  the  difference  being  more  noticeable 
for  harder  materials  (Xu,  2008). 

Another  investigation  into  the  effect  of  indenter  geometry  on  indentation  test  results 
was  performed  by  Sakharova  et  al.  (2009).  Conical  and  Berkovich  were  among  the  indenter 
geometries  examined.  Elastoplastic  indentations  using  Swift  hardening  law  (Wang,  2005) 
were  simulated  using  an  in-house  finite  element  code.  The  load-displacement  plots 
provided  in  the  paper  show  that  the  responses  from  the  different  indenters  were  very 
similar.  Small  differences  in  the  curves  were  observed  when  the  ratio  of  the  residual 
indentation  depth  after  reloading  and  the  indentation  depth  at  maximum  load  (hf/hmax) 
was  below  0.65.  In  these  cases  the  Berkovich  indenter  curves  were  slightly  higher 
(Sakharova,  2009). 

Comparisons  between  the  load-displacement  responses  of  materials  indented  by  conical 
and  Berkovich  indenters  using  finite  element  analysis  have  been  presented  in  these  earlier 
studies.  However,  the  results  do  not  seem  to  show  a  consistent  trend  and  therefore 
investigation  into  the  equivalence  between  these  two  types  of  indenters  is  still  incomplete. 
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II.  1 .3  Material  Modeling 


II.  1.3.1  Description  of  the  Constitutive  Model 

The  combined  isotropic/kinematic  hardening  model  is  a  pressure-independent  J2 
plasticity  model  that  can  be  used  for  modeling  materials.  The  behavior  of  the  material  is 
dependent  on  an  elastic  component  and  a  plastic  component.  To  show  this  the  strain  rate 
can  be  written  in  terms  of  a  plastic  and  elastic  component 

£  =  £el  +  £pl  (1.1) 
The  elastic  component  is  isotropic  and  linear.  It  is  dependent  of  two  material 
parameters  bulk  modulus,  K,  and  shear  modulus,  G.  The  parameters  are  calculated  based 
on  the  Young’s  modulus,  E,  and  Poisson’s  ratio, v 


K  = 


3(1  -2v) 


(1.2) 


G  = 


2(1  +  v) 


(1.3) 


Using  these  parameters  elasticity  can  be  decomposed  into  volumetric  and  deviatoric 
components 


1 

II 

a. 

(1.4) 

q=3Gsq 

(1.5) 

For  plastic  behavior  the  isotropic/kinematic  model  is  used  with  the  Mises  yield  surface 
in  ABAQUS  (Simulia,  2010).  The  yield  surface  is  given  by  the  function 


/  =  |(s-«)(s-«)-t2 


(1.6) 
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Figure  II.  1-1:  Yield  surface  of  combined  isotropic/kinematic  hardening  model. 


In  the  yield  surface  equation  S  is  the  deviatoric  stress  tensor  which  is  defined  as 
S  =  o  +  pi ,  a  is  the  deviatoric  component  of  the  backstress  tensor,  and  k  is  the  size  of  the 


yield  surface.  The  backstress  creates  a  translation  of  the  yield  surface  in  the  stress  space. 
This  shift  can  be  seen  in  the  preceding  figure  of  the  yield  surface.  A  Mises  yield  surface 
without  kinematic  hardening  would  remain  centered  along  the  hydrostatic  axis.  With  the 
kinematic  hardening  the  backstress  translates  the  yield  surface  so  that  the  center  is  no 
longer  coincident  with  the  hydrostatic  axis. 

Associated  plastic  flow  is  assumed  for  the  model  which  gives 


tPi  =  p  d/(g~«) 

da 


(1.7) 
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Where  is  the  equivalent  plastic  strain  rate  and  is  the  rate  of  plastic  flow. 

Hardening  consists  of  a  kinematic  component  and  an  isotropic  component.  The 
components  of  the  kinematic  hardening  are  the  backstresses,  of  which  there  can  be  several. 
The  hardening  law  for  each  backstress  is  given  as 


=  Ck  -j=  (a  -  a )s  pl  -  n aks  pl 


(1.9) 


N 

«=  X  a,  (1.10) 

k  =  1 

The  total  number  of  backstresses  is  N,  and  Ck  and  yk  are  material  parameters  that  are 
calibrated  in  ABAQUS  based  on  test  data.  The  values  Ck  are  the  initial  kinematic  hardening 
moduli,  and  yk  determine  the  rate  of  decrease  of  the  kinematic  hardening  moduli 
corresponding  to  an  increase  in  plastic  deformation  (Simulia,  2010). 

The  evolution  of  the  size  of  the  yield  surface,  k,  is  defined  by  the  isotropic  hardening 
behavior  of  the  model.  The  yield  surface  size  evolution  is  a  function  of  the  equivalent 

_ ryl 

plastic  strain  ,  s1,  and  can  be  written  as 

k  =  cx\n+QJ\-e-MPl)  (1.11) 


Where  cr  |0  is  the  size  of  the  yield  surface  when  plastic  strain  is  zero.  The  maximum 

change  in  the  size  of  the  yield  surface  is  given  by  Qx,  and  b  affects  the  rate  of  change  of  the 

yield  surface  size  while  plastic  strain  develops.  The  parameters  Qx  and  b  are  obtained 

through  calibration  with  test  data  (Simulia,  2010).  In  the  kinematic  hardening  model,  the 

kinematic  hardening  component  causes  a  shift  of  the  center  of  the  yield  surface  in  stress 
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space.  The  isotropic  component  can  expand  or  contract  the  yield  surface  range.  These 
components  allow  for  the  modeling  of  elastoplastic  behavior  of  materials. 

II.  1.3.2  Calibration  of  Material  Model 

The  experimental  results  and  numerical  simulations  presented  in  this  paper  were 
perfonned  on  6061  aluminum  as  a  typical  elastic-strain  hardening  plastic  material.  To 
evaluate  the  mechanical  properties  of  this  material,  a  macro  scale  uniaxial  test  was 
conducted  and  the  observed  stress-strain  curve  was  used  to  calibrate  the  constitutive  model 
used  in  the  finite  element  analysis.  All  analyses  were  conducted  using  ABAQUS  (Simulia, 
2010).  In  ABAQUS  a  metal  plasticity  constitutive  model  with  combined 
isotropic/kinematic  hardening  was  used  to  model  the  aluminum.  Calibration  of  the 
combined  isotropic/kinematic  hardening  model  in  ABAQUS  requires  half-cycle  stress- 
strain  test  data.  To  acquire  the  required  data  a  uniaxial  tension  test  was  performed  on  a 
“dog  bone”  specimen  of  6061  aluminum.  During  the  test  the  sample  was  loaded  in  uniaxial 
tension,  once  it  reached  initial  yield  stress  the  sample  was  unloaded  and  reloaded  to  yield 
again.  The  unloading  and  reloading  process  was  repeated  for  two  additional  cycles.  Stress 
and  strain  data  were  recorded  throughout  the  test. 

Table  II.  1-1:  Evolution  of  yield  stress  with  plastic  strain. 

Plastic  strain  0.0  0.002  0.01605  0.04667  0.07511 

Yield  stress  (kPa)  300,000  307,500  317,181  325,620  330,831 
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From  the  initial  elastic  loading  data  the  Young’s  modulus  (E)  was  calculated  as 
approximately  65GPa.  A  Poisson’s  ratio  (v)  of  0.33  was  assumed.  Given  the  value  of 
Young’s  modulus,  the  strain  hardening  data  relating  the  yield  stress  to  plastic  strain  were 
obtained  from  the  test  results.  Table  II.  1-1  shows  the  tabulated  values  of  yield  stress  versus 
plastic  strain.  This  data  was  used  to  calibrate  the  elastoplastic  constitutive  model. 

The  calibration  process  was  performed  by  conducting  a  plane  stress  simulation  of  a 
single  finite  element  subjected  to  displacements  that  approximate  those  applied  in  the 
uniaxial  tension  test.  After  conducting  a  series  of  simulations  and  adjusting  the  number  of 
back  stresses  to  five  for  smoothing  the  elastoplastic  transition,  the  finite  element  model 
provided  a  good  simulation  of  the  6061  aluminum  stress-strain  behavior,  as  shown  in 
Figure  II.  1-2. 
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Figure  II.  1-2:  Comparison  of  experimentally  obtained  stress-strain  curves  with  the 
simulation  of  a  combined  isotropic/kinematic  hardening  plasticity  model  for  6061 
aluminum. 


II.  1 .4  Finite  Element  Modeling 

In  this  study  three  different  finite  element  nanoindentation  models  were  investigated 
using  the  commercial  finite  element  analysis  program  ABAQUS.  The  first  model  was  a  2D 
axisymmetric  model  using  a  conical  indenter.  This  type  of  model  is  commonly  used  in 
research  and  allows  for  relatively  quick  computation.  The  second  model  was  a  3D  finite 
element  model  with  a  conical  indenter  tip.  This  model,  which  is  the  full  3D  equivalent  to 
the  2D  model,  was  used  to  confirm  the  suitability  of  the  3D  finite  element  mesh.  The  third 
model  was  a  3D  finite  element  mesh  that  included  complete  geometric  representation  of 
Berkovich  tip.  This  model  allows  for  a  thorough  comparison  of  Berkovich  tip  to  the 
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assumed  equivalent  conic  tip  advocated  in  literature.  Tip  geometries  for  all  the  above 
models  include  a  200nm  tip  radius  to  simulate  a  used  indenter  tip. 

II.  1 .4. 1  Verification  of  the  2D  Finite  Element  Model 

The  first  step  was  to  create  and  verify  a  2D  axisymmetric  nanoindentation  model  in 
ABAQUS.  The  perfonnance  of  the  developed  model  was  verified  against  the  available 
results  in  literature.  As  part  of  his  Ph.D.  thesis,  Poon  (2009)  perfonned  a  series  of  finite 
element  simulations  of  elastic  indentations  while  adjusting  the  radius  and  height  of  the 
specimen.  These  simulations  suggest  the  proper  specimen  size  to  approximate  a  body  of 
semi-infinite  extent  for  a  given  indentation  depth.  Reproducing  the  key  results  of  Poon’s 
simulation  helped  verify  the  quality  of  the  finite  element  mesh  and  the  sizing  of  the 
material  dimensions  used  in  the  present  study. 

Poon  (2009)  initially  provided  a  simulation  for  an  axisymmetric  indentation  onto  the 
top  surface  of  a  cylindrical  specimen  with  a  radius,  rs,  of  18pm  and  a  height,  hs,  of  30pm. 
The  specimen  was  modeled  as  an  elastic  isotropic  defonnable  solid  with  E=70GPa  and 
v=0.3.  The  70.3°  half-angle  “Berkovich  equivalent”  conical  indenter  was  modeled  as  an 
analytical  rigid  surface  with  a  tip  radius  of  p=200mn.  The  maximum  indentation  depth  was 
hmax=600nm.  The  finite  element  mesh  was  described  as  using  5006  three-node  linear 
axisymmetric  triangular  elements  (CAX3)  with  a  higher  mesh  density  at  the  indentation 
site.  For  later  simulations  the  number  of  nodes  was  scaled  depending  on  the  specimen 
dimensions.  Contact  between  the  indenter  and  specimen  was  modeled  as  frictionless. 
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Figure  II.  1-3:  Face  partitioning  of  2D  axisymmetric  specimen. 


A  2D  axisymmetric  model  was  created  for  the  present  study  using  the  same  parameters 
detailed  in  Poon  (2009).  While  CAX3  element  type  was  used,  a  different  finite  element 
mesh  was  constructed  to  see  the  potential  effects  of  mesh  configuration.  The  specimen  was 
modeled  with  editable  dimension  constraints  for  the  radius  and  height.  The  face  of  the 
specimen  was  then  partitioned  into  quarters.  In  order  to  easily  provide  a  higher  density 
mesh  at  the  indentation  site  the  quadrant  containing  the  indentation  site  (upper  left)  was 
quartered.  This  last  step  was  repeated  until  the  face  in  the  comer  of  the  indentation  site  was 
sufficiently  small.  The  resulting  face  partition  of  the  specimen  is  shown  in  Figure  II.  1-3. 
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Edges  were  then  seeded  to  provide  a  high  mesh  density  at  the  indentation  site  and  less 
density  further  away  as  shown  in  Figure  II.  1-4. 


Figure  II.  1-4:  2D  axisymmetric  finite  element  mesh  of  the  specimen. 


A  linear  elastic  analysis  was  perfonned  using  the  same  material  size  and  specifications 
used  in  the  initial  Poon  simulation.  Figure  II.  1-5  shows  the  load-displacement  result 
obtained  from  this  simulation  and  compares  it  to  the  simulation  result  presented  in  Poon 
(2009).  The  curves  were  in  good  agreement  with  only  slight  differences  which  may  be 
caused  by  errors  in  reading  the  values  from  the  curves  in  Poon  (2009). 


33 


Figure  II.  1-5:  Load-displacement  comparison  for  elastic  indentation  simulation  into  a 
specimen  with  rs=18pm  and  hs=  30pm. 


More  verification  simulations  were  conducted  by  varying  the  tip 
The  specimen  in  these  simulations  was  described  as  having 
Converged  geometry  is  determined  by  using  the  following  relations: 


h 

— ^  >100 


radius  of  the  indenter. 
converged  geometry. 


(1.12) 

(1.13) 
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The  result  is  shown  for  the  range  of  indentation  of  approximately  260mn,  though  the  actual 
indentation  depth  is  not  given.  To  reproduce  the  results  for  this  study  an  indentation  depth 
of  300nm  was  used  on  a  material  with  rs=  30pm  and  hs=  30pm.  Load-displacement  results 
are  shown  in  Figure  II.  1-6.  The  results  show  good  agreement  in  values  and  trends  when 
compared  with  those  obtained  by  Poon  (2009). 
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Figure  II.  1-6:  Load-displacement  comparison  for  elastic  indentation  simulation  into  a  specimen 
with  varying  indenter  tip  radii,  p  (30,  75,  120,  150  and  200  nm). 


Simulations  were  also  performed  keeping  the  same  specimen  radius  of  rs=18pm,  while 
varying  the  specimen  height,  hs,  through  the  values  18,  30,  40,  60,  92,  and  120pm.  The 
load-displacement  results  for  the  simulations  are  shown  in  Figure  II.  1-7.  The  results  exhibit 
the  expected  trend  where  load  increases  with  decreased  specimen  height. 
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For  further  simulations  in  the  study  the  material  dimensions  are  chosen  to  be  rs=  60pm  and 
hs=  60pm  to  accommodate  indentations  of  up  to  600nm  according  the  convergence 
relations  provided. 


Figure  II.  1-7:  Load-displacement  curves  for  elastic  indentation  simulation  into  a  specimen  with 
rs  =  18  pm  and  varying  hs  (18,  30,  42,  60,  92  and  120  pm). 


II.  1.4.2  Two  and  Three  Dimensional  Elastoplastic  Simulations 

The  2D  axisymmetric  finite  element  mesh  was  verified  for  elastic  indentation.  To 

perform  elastoplastic  indentations  the  calibrated  elastoplastic  constitutive  model  was  used. 

The  material  specimen  was  sized  to  accommodate  indentations  up  to  a  600nm.  The  element 

type  in  the  mesh  was  changed  to  a  four-node  bilinear  axisymmetric  quadrilateral  element 

36 


(CAX4R),  which  for  linear  elastic  analysis  would  produce  the  same  results  as  those 
obtained  by  the  CAX3  element  used  in  verification  studies  discussed  earlier. 

Initial  elastoplastic  simulations  resulted  in  load-displacement  curves  that  exhibited  a 
slight  oscillation  or  “stepping  type”  behavior  in  the  load-displacement  curve.  This  behavior 
was  corrected  by  providing  a  more  unifonn  grid  for  the  fine  mesh  around  the  indentation 
site.  The  mesh  was  verified  through  a  convergence  study  of  the  2D  elastoplastic  indentation 
as  shown  in  Figure  II.  1-8.  The  final  mesh  chosen  included  8,582  elements.  The  2D  analysis 
results  were  later  used  as  a  baseline  for  comparison  against  the  3D  indentation  models. 
Figure  II.  1-4  shows  the  2D  axisymmetric  mesh  with  a  detail  of  the  local  mesh  near  the 
indenter  tip. 


Figure  II.  1-8:  Comparison  of  the  simulations  using  different  2D  axisymmetric  finite 
element  meshes  in  elastoplastic  indentation  for  606 1  aluminum. 
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To  provide  confirmation  of  the  suitability  of  3D  meshes  and  for  direct  comparison  to 
the  2D  model,  a  3D  conical  tip  nanoindentation  model  was  first  created.  Similar  to  the  2D 
model  the  3D  indentation  model  was  created  using  a  70.3°  half-angle  Berkovich  equivalent 
cone.  The  conical  indenter  was  created  as  a  revolved  shell  analytical  rigid  surface  as  shown 
in  Figure  II.  1-9. 


Figure  II.  1-9:  3D  model  of  70.3°  half-angle  conical  indenter. 


The  specimen  was  created  in  a  cylindrical  shape  with  the  same  depth  and  radius  as  the 
2D  axisymmetric  model  and  elastoplastic  material  parameters  were  also  identical  to  those 
used  in  the  2D  axisymmetric  model.  The  finite  element  mesh  of  the  specimen  was  created 
so  that  it  was  fine  on  the  surface  around  the  site  of  indentation  as  well  as  the  initial  depth  of 
the  cylinder.  The  mesh  was  then  increasingly  coarsened  in  the  areas  that  were  further  away 
from  the  indentation  site.  Each  quadrant  of  the  specimen’s  top  surface  was  meshed  by 
using  a  method  similar  to  that  used  in  the  2D  axisymmetric  model,  i.e.,  creating 
increasingly  smaller  square  partitions  and  then  seeding  the  edges  for  desired  mesh  density. 
The  element  type  was  chosen  to  be  an  8-node  linear  brick  using  reduced  integration  and 
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hourglass  control  (C3D8R).  Figure  II.  1-10  shows  the  3D  finite  element  mesh  and  a  local 
mesh  detail  near  the  tip  is  shown  in  Figure  II.  1-1 1 . 


Figure  II.  1-10:  3D  finite  element  mesh  of  the  specimen. 


Figure  II.  1-1 1:  Local  details  of  3D  finite  element  mesh  near  the  indenter. 
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A  3D  Berkovich  tip  nanoindentation  model  was  the  third  and  final  finite  element  mesh 
created  for  analysis  and  comparison.  The  material  model  and  the  finite  element  mesh  of  the 
specimen  were  carried  over  from  the  conical  simulation.  The  main  difference  between  the 
two  3D  models  was  the  tip  geometry.  The  Berkovich  tip  was  modeled  as  a  three  sided 
pyramid  with  a  semi-angle  between  vertical  and  each  face  of  65.27°.  An  area  of  difficulty 
in  modeling  the  Berkovich  tip  was  in  the  process  of  rounding  the  tip.  A  pyramidal  shape 
cannot  be  tangentially  rounded  by  a  sphere  along  all  sides  and  edges.  Instead  the  tip  was 
rounded  using  the  intersection  between  a  sharp  Berkovich  geometry  and  a  conical 
geometry  with  a  200nm  tip  radius  and  a  half  angle  of  77.05°  which  is  equivalent  to  the 
semi-angle  to  the  sharp  edges  of  a  Berkovich  pyramid.  Due  to  the  complexity  in  the 
Berkovich  tip  geometry  it  cannot  be  treated  as  an  analytical  rigid  surface,  instead  it  was 
modeled  as  a  discrete  rigid  surface.  A  3D  model  of  the  Berkovich  indenter  tip  is  shown  in 
Figure  II.  1-12. 


Figure  II.  1-12:  3D  model  of  Berkovich  indenter. 


After  some  preliminary  simulations,  slight  adjustments  were  made  to  the  specimen 


mesh  for  the  Berkovich  indentation  to  accommodate  interaction  with  the  more  complicated 
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Berkovich  tip  geometry.  Distortions  in  the  mesh  would  arise  depending  on  the  orientation 
of  elements  that  were  contacted  by  the  pyramid  edge.  The  mesh  was  modified  in  the  region 
of  tip  contact  to  align  element  edges  with  the  pyramid  edges.  This  correction  reduced  mesh 
distortion  and  provided  an  improvement  in  the  load-displacement  response.  A  local  detail 
of  the  improved  mesh  configuration  is  shown  in  Figure  II.  1-13. 


Figure  II.  1-13:  Local  top  surface  details  of  the  modified  3D  finite  element  mesh 
configured  to  improve  interaction  with  Berkovich  indenter. 


All  three  finite  element  models  shared  several  common  characteristics.  The 
elastoplastic  material  properties  used  were  those  obtained  from  the  uniaxial  tension  test  of 
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6061  aluminum.  Contact  was  modeled  as  frictionless  between  the  indenter  and  material. 


Displacement  of  the  indenter  was  applied  through  a  displacement  boundary  condition 
applied  to  a  reference  node  of  the  rigid  indenter  surface.  Indentation  force  was  obtained  by 
requesting  output  of  the  reaction  force  at  the  reference  node  during  the  analysis. 

Results  in  the  study  were  all  obtained  from  ABAQUS  implicit  analysis.  However, 
explicit  analysis  was  also  considered.  Simulation  results  revealed  that  explicit  analysis 
would  work  well  for  analysis  of  a  3D  cone  since  it  provided  a  response  very  similar  to  that 
of  the  implicit  analysis.  However,  the  explicit  analysis  provided  poor  results  for  the 
Berkovich  indentation.  The  poor  performance  was  attributed  to  large  elemental  distortion 
of  the  material  in  regions  of  the  mesh  that  had  contact  interaction  with  the  edges  of  the 
Berkovich  pyramid.  This  led  to  a  decrease  in  simulation  accuracy  with  increased  mesh 
refinement  due  to  larger  relative  distortions  for  the  finer  elements.  These  large  distortions 
not  only  led  to  inaccuracy  but  the  unexpected  behavior  where  the  material  would  pull  back 
on  the  indenter  resulting  in  a  negative  reaction  force  during  unloading.  Figure  II.  1-14 
shows  an  example  of  the  mesh  distortion  experienced  in  an  explicit  analysis.  Implicit 
analysis  was  confirmed  to  be  the  preferable  method  for  nanoindentation  simulation. 
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Figure  II. 1-14:  Mesh  distortion  observed  in  explicit  simulation  of  Berkovich  indentation. 


Using  these  three  finite  element  models,  a  comparison  of  how  geometry  affects  the 
indentation  results  can  be  made.  The  main  point  of  interest  was  how  similar  the  load- 
displacement  behavior  was  between  a  Berkovich  indenter  and  the  equivalent  conical 
indenter. 
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II.  1.5  Results  and  Discussion 


II.  1.5.1  Finite  Element  Simulations  of  Nanoindentation 

The  results  were  first  examined  by  considering  the  load-displacement  curves  of  the 
elastoplastic  indentation  process.  The  load-displacement  curve  represents  the  overall 
material  reaction  to  the  nanoindentation  process  and  it  also  represents  what  would  be 
recorded  during  a  nanoindentation  experiment.  An  indentation  depth  of  400mn  was  chosen 
in  this  comparison  study. 

Figure  II.  1-15  compares  the  load-displacement  results  obtained  from  a  2D 
axisymmetric  conical  indentation  simulation  and  a  3D  conical  indentation  simulation.  It 
would  be  expected  that  both  simulations  would  provide  the  same  curve.  In  line  with 
expectations  the  two  curves  are  in  good  agreement  with  each  other.  At  maximum 
indentation  depth  the  difference  between  the  two  curves  is  less  than  2%.  The  differences 
may  be  attributed  to  the  fact  that  the  mesh  of  the  2D  axisymmetric  model  is  much  finer 
than  the  3D  conical  model  in  the  area  around  the  tip.  Increasing  the  mesh  refinement  in  the 
3D  model  would  lead  to  significant  increases  in  simulation  time  while  producing  little 
improvement  in  the  simulation  results.  Hence,  the  3D  finite  element  mesh  was  considered 
to  be  sufficiently  refined. 

Figure  II.  1-16  shows  a  comparison  of  the  load-displacement  curves  for  the  Berkovich 
and  conical  indenters.  It  can  be  seen  that  there  is  a  difference  between  the  results.  The 
Berkovich  indentation  curve  is  noticeably  lower  than  that  of  the  conical  indenter. 
According  to  common  practice  the  curves  should  match  if  the  projected  area  to  indentation 


depth  is  equal.  The  contact  area  during  indentation  was  then  checked  for  both  indenter 
geometries  to  determine  if  there  are  differences  in  contact  area  that  are  leading  to  load 
differences.  Figure  II.  1-1 7  shows  the  contact  area  for  both  Berkovich  and  conical  indenters 
as  it  relates  to  indentation  depth  in  the  simulations.  The  contact  area  curves  for  the  two 
indenters  have  slight  differences  but  in  general  they  are  very  similar  to  each  other. 


Figure  II.  1-15:  Comparison  of  load-displacement  curves  for  2D  and  3D  conical 
indentation  simulations  in  606 1  aluminum. 
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Figure  II.  1-16:  Comparison  of  load-displacement  curves  for  conical  and  Berkovich 
indentation  simulations  in  606 1  aluminum. 
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Figure  II.  1-17:  Comparison  of  contact  areas  computed  in  conical  and  Berkovich 
indentation  simulations  in  606 1  aluminum. 

The  contact  areas  for  Berkovich  and  conical  indenters  are  essentially  the  same  but  there 
load-displacement  curves  have  noticeable  differences;  this  provides  good  motivation  to 
examine  the  normalized  contact  stresses  for  the  two  indenters.  The  normalized  contact 
stresses  during  indentation  for  the  two  indenters  is  shown  in  Figure  II.  1-18.  Both  curves 
share  the  same  initial  elastic  loading  during  contact  initiation;  however  after  the  transition 
into  the  plastic  regime  the  curve  for  the  normalized  stress  for  the  Berkovich  indenter  is 
consistently  lower  than  the  curve  for  the  conical  indenter.  There  seems  to  be  a  geometry 
effect  that  is  resulting  in  different  indentation  responses  for  the  two  indenters.  It  is  also 
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observed  that  the  contact  stress  is  approximately  2.5  times  the  yield  stress.  This  is  possible 
because  the  constitutive  model  is  shear  sensitive  and  not  pressure  sensitive. 


Figure  II.  1-1 8:  Comparison  of  nonnalized  contact  stresses  computed  for  conical  and 
Berkovich  indentation  simulations  in  6061  aluminum. 


To  further  examine  the  potential  effect  of  geometry  during  the  elastoplastic  indentation, 

finite  element  visualizations  were  utilized  to  examine  the  material  behavior  around  the 

indentation  site.  Figure  II.  l-19(a)  shows  the  Mises  stress  field  for  the  70.3°  conical  indenter 

at  maximum  indentation  depth.  The  axisymmetric  geometry  of  the  conical  indenter 

provides  an  axisymmetric  stress  field.  Figure  II.  1-  19(b)  shows  the  Mises  stress  field  for  the 

Berkovich  indenter  with  a  view  cut  in  the  plane  perpendicular  to  a  sharp  edge  of  the 
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pyramid  which  provides  an  angle  similar  to  that  of  the  equivalent  cone.  On  this  cross- 
section  the  stress  fields  between  the  two  indenters  are  quite  similar.  Figure  II.  1- 19(c)  shows 
the  Mises  stress  field  for  the  Berkovich  indenter  with  a  view  cut  made  in  the  plane  of  a 
sharp  edge  of  the  pyramid  which  provides  the  cross  section  with  the  most  dissimilar  angle. 
On  this  cross-section  the  stress  field  extends  noticeably  further  laterally  on  the  side  of  the 
flat  edge  of  the  pyramid  and  does  not  extends  as  far  near  the  sharper  edge. 

The  distribution  of  equivalent  plastic  strain  which  depends  on  the  material  model  was 
also  examined.  Figure  II.l-20(a)  shows  the  equivalent  plastic  strain  distribution  on  the 
material  surface  after  conical  indentation.  The  maximum  equivalent  plastic  strain  for  the 
conical  indentation  was  170%.  The  distribution  of  equivalent  plastic  strain  for  the 
Berkovich  indentation  is  shown  in  Figure  II.  1 -20(b).  For  the  Berkovich  indentation  the 
maximum  equivalent  plastic  strain  was  204%  for  the  same  indentation  depth  as  that  in 
conical  indentation.  It  can  be  seen  from  the  strain  distribution  that  the  edges  of  the  pyramid 
contribute  to  creating  higher  equivalent  plastic  strains  in  the  Berkovich  indenter  than  those 
experienced  by  the  conical  indenter. 
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Figure  II.  1-19:  Mises  stress  in  6061  aluminum  for  (a)  70.3°  conical  indenter,  (b)  Berkovich 
indenter,  viewed  in  plane  perpendicular  to  indenter  edge,  and  (c)  Berkovich  indenter, 
viewed  in  plane  of  indenter  edge. 
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Figure  II. 1-20:  Contours  of  equivalent  plastic  strain  in  6061  aluminum  for  (a)  70.3°  conical 
indenter  and  (b)  Berkovich  indenter. 
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Despite  the  contact  area  being  the  same  throughout  the  indentation  process  for  both 
indenters,  the  material  is  responding  differently  due  to  geometry  effects.  The  edges  of  the 
Berkovich  pyramidal  geometry  induce  higher  plastic  strains.  This  means  that  for  the 
Berkovich  indenter  there  is  more  localized  yielding  in  the  material.  With  more  localized 
yielding  the  material  is  less  resistant  to  penetration  which  results  in  lower  required 
indentation  loads  for  the  Berkovich  indenter  compared  to  the  70.3°  conical  indenter. 

II.  1.5.2  Comparison  of  Nano  indentation  Simulations  and  Experiments 

To  assess  the  validity  of  numerical  simulations,  it  is  always  desirable  to  compare  them 
with  the  observed  responses  of  the  material  in  nanoindentation  experiments.  This  will  also 
provide  additional  confidence  regarding  the  ability  of  the  finite  element  analysis  to 
accurately  predict  a  real  nanoindentation  test  result.  Hence,  a  cylindrical  specimen  was 
prepared  for  the  nanoindentation  test  using  the  same  type  of  aluminum  used  in  calibration 
of  the  finite  element  model.  A  series  of  nanoindentation  test  were  perfonned  using  an 
Asylum  Research  MFP-3D™  Stand  Alone  Atomic  Force  Microscope  fitted  with  a 
Nanolndenter  module  (Asylum  Research,  Santa  Barbara,  CA).  A  series  of  64  indentations 
were  performed  in  a  grid  pattern  on  the  material  surface  and  a  sampling  of  curves  showing 
the  variation  in  load-displacement  response  among  the  indentations  is  shown  in  Figure  II.  1- 
2 1 .  The  indentations  were  perfonned  using  force  control  to  a  maximum  indentation  force 
of  1.4mN.  The  variation  in  indentation  depth  is  thought  to  be  from  variations  in  the 
material  properties  due  to  the  effects  of  machining  the  sample.  An  AFM  produced  contour 
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plot  of  the  material  surface  is  shown  in  Figure  II.  1-22.  The  surface  roughness  for  this 
100pm2  plot  was  given  as  approximately  3.5nm. 


Figure  II.  1-2 1 :  Upper  and  lower  bounds  of  experimental  load-displacement  curves  in  6061 
aluminum. 


Figure  II.  1-22:  AFM-produced  contour  plot  of  6061  aluminum  specimen  surface. 
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Finite  element  simulations  were  performed  to  a  maximum  depth  similar  to  the 
experimental  indentations.  Comparing  the  unloading  slopes  of  the  experimental  and 
simulated  curves  revealed  a  Young’s  modulus  from  the  nano  indentation  experiments  of 
approximately  59GPa,  which  was  approximately  9%  lower  in  magnitude  than  the  65GPa 
modulus  detennined  in  the  macro  scale  test.  The  simulation  results  for  both  the  conical  and 
Berkovich  indenters  provided  load-displacement  curves  with  significantly  lower  loads 
compared  to  experimental  values.  Warren  et  al.  (2006)  discuss  that  machining  induced 
residual  stress  can  significantly  affect  load-displacement  curves.  The  residual  stresses  may 
be  estimated  through  finite  element  simulation.  In  their  sample  of  52100  steel,  residual 
stresses  were  estimated  to  be  18%  of  yield  stress  in  the  horizontal  direction  and  12%  of 
yield  stress  in  the  vertical  direction.  These  values  were  used  as  a  starting  point  in  estimating 
the  residual  stresses  in  the  aluminum  specimen  in  the  present  study. 

Using  the  2D  axisymmetric  finite  element  model  (due  to  its  lower  computational  effort) 
a  series  of  simulations  were  perfonned  while  adjusting  residual  stress  parameters.  The 
residual  stresses  were  estimated  to  be  approximately  8%  of  yield  stress  in  the  horizontal 
direction  and  5%  of  yield  in  the  vertical  direction.  Incorporating  these  residual  stresses  as 
the  initial  stresses  was  shown  to  bring  up  the  computed  load-displacement  curves  well 
within  the  experimental  bounds.  Figure  II.  1-23  compares  the  conical  and  Berkovich 
indentations  without  residual  stresses,  the  conical  indentation  with  residual  stresses,  and  the 
experimental  bounds. 
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Figure  II.  1-23:  Comparison  of  numerically  obtained  force-displacement  curves  with 
experimental  bounds  for  606 1  aluminum. 


II.  1 .6  Conclusions 

Finite  element  simulations  of  elastoplastic  nano  indentations  were  performed  for  conical 
(70.3°)  and  Berkovich  indenters  to  study  the  effect  of  indenter  geometry  on  the  response  of 
an  elastoplastic  material.  It  is  shown  that  there  are  clear  differences  in  the  load- 
displacement  responses  between  the  two  indenters.  Loads  in  the  Berkovich  indentation  are 
noticeably  lower.  It  is  confirmed  that  contact  areas  during  indentation  are  very  similar, 
resulting  in  lower  normalized  contact  stress  for  the  Berkovich  indenter.  Examination  of 
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stress  fields  near  the  indenter  tip  shows  that  the  symmetrical  cross-sections  of  the 
Berkovich  indenter  are  similar  to  the  axisymmetric  stress  field  of  the  conical  indenter  and 
differences  arise  when  compared  to  asymmetric  cross-sections  of  the  Berkovich  indenter. 
Maximum  equivalent  plastic  strain  induced  by  indenter  is  demonstrated  to  be  20%  higher 
for  the  Berkovich  indenter.  It  is  concluded  that  higher  localized  plastic  strains  would  result 
in  less  resistance  to  indentation  by  the  material  and  explain  the  load-displacement  behavior. 
The  results  reveal  that  the  use  of  a  70.3°  conical  indenter  to  approximate  the  response  from 
a  Berkovich  indenter  is  not  valid  for  all  materials. 

The  ability  of  finite  element  analysis  to  predict  experimental  nanoindentation  results 
was  also  examined.  Comparisons  with  experimental  results  suggest  that  machining-induced 
residual  stresses  have  likely  affected  the  experimental  results.  It  is  also  demonstrated  that 
these  residual  stresses  may  be  estimated  by  matching  simulated  load-displacement  curves 
within  the  range  of  experimental  curves. 
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II. 2:  Nanoindentation  in  a  Viscoelastic  Material 


II.2.1  Introduction 

Nanoindentation  in  viscoelastic  materials  presents  new  challenges  due  to  the  time 
dependent  behavior  of  the  material.  Thus,  the  response  of  the  material  is  sensitive  to  the 
speed  of  penetration  by  the  indenter.  With  elastic  and  elastoplastic  materials  it  is  common 
for  material  parameters  such  as  the  elastic  modulus  to  be  detennined  by  examining  the 
unloading  curve  obtained  during  a  nanoindentation  test  (Oliver,  2004).  The  techniques  are 
well  established  and  work  well  for  elastic  and  elastoplastic  materials.  However,  the 
procedures  developed  for  these  materials  are  either  inaccurate  or  have  many  uncertainties 
for  viscoelastic  materials  (VanLandingham,  2001).  The  time  dependence  adds  complexities 
to  the  material  behavior  that  the  previous  techniques  cannot  capture.  For  instance,  the 
Young’s  modulus  of  a  material  is  no  longer  constant  like  it  is  with  elastic  and  majority  of 
elastoplastic  materials. 

There  is  still  much  research  needed  in  the  investigation  of  viscoelastic  nanoindentation. 
The  use  of  finite  element  analysis  provides  the  option  to  quickly  examine  the  behavior  of 
this  type  of  material  by  varying  key  parameters  and  observing  the  differences  in  behavior. 
The  numerical  simulations  also  provide  the  ability  to  look  more  closely  at  the  behavior  of  a 
material  in  nanoindentation.  However,  for  finite  element  simulations  to  be  useful  they  must 
also  be  proven  to  accurately  predict  the  behavior  observed  in  experimental  nanoindentation 
curves. 

In  this  study  the  goal  is  to  examine  how  well  a  viscoelastic  nanoindentation  simulation 
can  capture  the  behavior  of  the  time  dependent  material.  First  it  will  be  confirmed  that  a 


viscoelastic  material  model  can  be  calibrated  given  experimental  relaxation  data.  This 
material  response  can  be  defined  using  Prony  series  parameters.  Then,  the  study  will 
investigate  the  nanoindentation  response  of  the  viscoelastic  material  and  its  comparison  to 
experimental  nanoindentation  results  for  the  same  type  of  material. 

II.2.2  Literature  Review 

An  examination  of  depth-sensing  indentation,  as  used  to  measure  the  elastic  modulus  of 
polymeric  materials,  was  performed  by  VanLandingham  et  al.  (2001).  In  their  paper  the 
authors  reviewed  the  traditional  analysis  of  load-displacement  data  for  nanoindentation  and 
discuss  the  limitations  and  uncertainties  that  are  found  for  various  measurements  in 
nanoscale  indentation.  The  authors  also  discuss  that  the  application  of  these  established 
methods  toward  use  in  polymers  often  lead  to  the  inaccurate  measurements  of  the  elastic 
modulus.  The  paper  shows  that  in  comparison  to  bulk  measurements  obtained  in 
macroscopic  experiments  the  elastic  modulus  values  measured  in  quasi-static  indentations 
are  relatively  high.  The  differences  are  attributed  to  the  effect  of  the  viscoelastic  behavior 
in  the  material  and  its  effect  on  the  shape  of  the  load-displacement  unloading  curve.  As  a 
method  to  potentially  alleviate  many  of  the  problems  associated  with  the  quasi-static 
indentation  testing  VanLandingham  et  al.  (2001)  discuss  the  use  of  dynamic  indentation 
testing.  The  authors  discuss  that,  previously,  rigorous  analysis  of  dynamic  indentation 
behavior  of  polymers  with  regard  to  viscoelasticity  has  not  been  reported.  The  study 
concludes  that  for  advancements  in  quantitative  characterization  of  polymer  properties 
there  will  need  to  be  material-independent  calibration  procedures,  reference  materials, 
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advances  in  instrumentation,  and  new  analysis  and  testing  procedures  to  account  for 
polymer  viscoelastic  behavior  (VanLandingham,  2001). 

A  study  by  Cheng  and  Cheng  (2005)  looks  at  conical  indentation  in  linear  viscoelastic 
materials  and  the  relationship  between  contact  depth,  initial  unloading  slope,  and  the 
viscoelastic  properties  of  the  material.  In  particular  they  examine  the  Oliver-Pharr  method 
for  determining  the  contact  depth  and  whether  this  method  is  applicable  to  indentations  in 
viscoelastic  materials.  In  their  study  the  authors  use  finite  element  analysis  to  examine  the 
viscoelastic  behavior.  The  study  detennines  a  relationship  between  initial  unloading  slope, 
contact  depth,  and  the  instantaneous  modulus  for  a  sufficiently  high  rate  of  unloading.  A 
relationship  between  contact  depth  and  indenter  displacement  is  also  discussed.  From  their 
results  the  authors  conclude  that  the  Oliver-Pharr  method  of  estimating  contact  depth  is  not 
applicable  for  nanoindentation  of  viscoelastic  materials.  These  conclusions  are  based  on 
constant  displacement  rate  loading  and  are  said  to  hold  for  other  loading  conditions  (Cheng 
and  Cheng,  2005). 

Other  studies  have  also  looked  into  the  accuracy  of  the  Oliver-Pharr  method  to 
determine  material  properties  of  viscoelastic  materials.  Tranchida  et  al.  (2006)  looked  into 
the  Oliver-Pharr  method’s  ability  to  accurately  evaluate  Young’s  modulus  of  polymers  by 
nanoindentation.  When  comparing  Young’s  modulus  detennined  in  macroscopic  tests  to 
the  modulus  determined  through  the  Oliver-Pharr  method  for  a  viscoelastic  material  the 
study  shows  that  there  is  clear  disagreement  between  the  values.  Young’s  modulus  is 
shown  to  be  overestimated  up  to  3.3  times  by  the  Oliver-Pharr  method.  The  study  proposes 
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the  use  of  a  correction  factor  relating  the  modulus  obtained  through  nanoindentation  and 
the  expected  value.  The  authors  conclude  that  a  satisfactory  evaluation  of  the  Young’s 
modulus  for  polymers  can  be  obtained  by  coupling  their  correction  factor  with  the 
calibrations  of  the  Oliver-Pharr  method  (Tranchida,  2006). 

In  another  study  looking  further  into  the  issue  of  using  the  Oliver-Pharr  method  for 
detennining  Young’s  modulus,  Tranchida  et  al.  (2007)  examines  the  effect  of  pile-up  and 
viscoelasticity  on  the  overestimation  of  Young’s  modulus.  The  study  shows  that  the  effect 
of  pile-up  can  be  minimized  by  using  shallow  indentations  at  a  true  nanoscale.  However, 
evaluating  the  Young’s  modulus  with  the  Oliver-Pharr  method  still  gives  rise  to  significant 
deviations  with  respect  to  the  value  measured  macroscopically.  One  issue  with  viscoelastic 
indentations  is  the  emergence  of  a  “nose”  in  the  unloading  force  curve.  This  can  be 
minimized  by  perfonning  indentation  at  large  loading  rates.  Again  the  Young’s  modulus 
detennined  through  the  Oliver-Pharr  method  is  still  concluded  to  be  consistently  higher 
than  the  bulk  elastic  modulus.  It  appears  that  viscoelasticity  has  a  significant  effect  on  the 
contact  mechanics  compared  to  the  elastic  one.  The  conclusion  is  that  the  Oliver-Pharr 
procedure  and  other  methods  derived  from  elastic  contact  cannot  be  applied  to  the  load- 
displacement  curve  of  viscoelastic  nanoindentations  (Tranchida,  2007). 

Nanoindentation  is  often  used  in  the  detennining  of  the  material  properties  in  thin 
films.  Zhou  and  Komvopoulos  (2006)  examine  the  viscoelastic  behavior  of  thin  films, 
showing  that,  with  thin  films,  stiffness  increases  with  the  indentation  depth  due  to  increased 
contact  area  and  the  effect  of  the  rigid  substrate.  Thickness  and  loading  rate  are  discussed 
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as  having  a  strong  effect  on  the  films  stiffness.  Similar  to  other  studies  the  Oliver-Pharr 
method  of  detennining  Young’s  modulus  overestimates  its  value  though  with  the  thin  films 
the  substrate  plays  a  role  in  the  overestimation.  The  study  shows  that  experimental  results 
for  reduced  elastic  modulus  which  were  modified  by  subtracting  the  viscous  component 
show  fair  agreement  with  analytical  results  (Zhou  and  Komvopoulos,  2006). 

Spherical  indentation  of  linear  viscoelastic  materials  are  examined  in  a  paper  by 
Ramesh  Kumar  and  Narasimhan  (2004).  In  this  study  the  relaxation  response  of 
polymethyl  methacrylate  (PMMA)  detennined  from  a  conventional  mechanical  test.  This 
data  from  the  relaxation  response  can  be  used  to  detennine  the  Prony  series  which  is  an 
exponential  series  that  defines  the  relaxation  function  of  the  material.  The  study  assesses 
the  validity  of  viscoelastic  indentation  theories  by  comparing  experimental  and  finite 
element  indentation  results.  For  spherical  indentations  the  two  are  shown  to  be  quite 
similar.  The  authors  show  that  depth-sensing  indentation  may  be  used  to  extract 
viscoelastic  response  function  but  that  the  conventional  method  of  using  the  unloading 
curve  may  results  in  large  errors  in  the  detennination  of  Young’s  modulus  (Ramesh 
Kumar,  2004). 

Another  study  that  includes  examination  of  the  indentation  response  of  polymethyl 
methacrylate  (PMMA)  was  perfonned  by  Oyen  (2006).  In  the  study  indentations  were 
perfonned  on  PMMA  using  a  Berkovich  indenter.  It  is  shown  that  additional  displacements 
occurred  beyond  that  due  to  the  viscoelastic  response.  These  additional  displacements  are 
likely  due  to  plastic  defonnations  and  are  mainly  caused  by  the  sharp  edges  of  the 
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Berkovich  indenter  and  complicate  the  viscoelastic  analysis.  The  author  states  that  a 
spherical  indenter  may  be  a  better  choice  to  keep  the  material  response  within  the  linear 
viscoelastic  regime.  If  pyramidal  indenter  geometry  is  used  then  plastic  defonnations  must 
be  considered  (Oyen,  2006). 

Sane  and  Knauss  (2001)  perfonned  a  study  to  detennine  various  material  functions  of 
PMMA.  The  paper  presents  data  for  the  relaxation  response  of  PMMA  which  are  of 
particular  interest  to  this  thesis.  This  provides  an  opportunity  to  model  a  second  PMMA 
material,  in  addition  to  that  reported  by  Ramesh  Kumar  and  Narasimhan  (2004),  and 
provide  an  opportunity  to  make  a  comparison  between  the  indentation  responses  of  two 
different  reported  samples  of  PMMA. 

There  have  been  many  studies  on  the  response  of  viscoelastic  materials  and  their 
response  during  indentation.  There  is  a  clear  consensus  that  conventional  method  used  to 
detennine  properties  of  elastic  and  elastoplastic  materials  are  not  accurate  for  viscoelastic 
materials.  Further  investigation  and  insights  are  needed  in  regards  to  the  response  of  these 
materials  during  nanoindentation. 

II.2.3  Material  Modeling 

II.2.3.1  Description  of  the  Constitutive  Model 

The  viscoelasticity  model  describes  rate-dependent  material  behavior.  These  materials 
experience  dissipative  losses  that  are  caused  by  internal  damping  effects.  Considering  small 
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strains,  let  us  assume  that  that  material  is  subjected  to  a  time  varying  shear  strain,  y(t).  The 
shear  stress  response,  x(t),  is  defined  as 

T(t)  =  ^GR(t-s)y{s)ds  (2.1) 

Where  the  time  dependent  shear  relaxation  modulus,  given  as  GR(t),  characterizes  the 
material  response.  The  constitutive  behavior  is  clearly  demonstrated  by  a  relaxation  test  in 
which  a  sudden  strain  y  is  applied  to  the  material  specimen  and  is  then  held  constant.  At 
zero  time,  taken  as  when  the  strain  is  suddenly  applied,  the  shear  stress  is  given  as 

r(t)  =  GR(t)r  (2.2) 

Since  y  is  held  constant  as  a  fixed  strain,  y  =0  for  all  t>0,  after  a  long  period  of  time  the 

material  will  converge  to  a  constant  stress  (Simulia,  2010). 

In  ABAQUS  the  material  relaxation  is  defined  using  the  dimensionless  form  given 
below: 

gR(t)  =  GR(t)/  G0  (2.3) 

The  instantaneous  shear  modulus  is  G0  =  G^(0).  Given  the  dimensionless  form  the 
equation  for  shear  stress  becomes: 

T(t)  =  G0\[gR{t-s)y{s)ds  (2.4) 

The  dimensionless  relaxation  function  gR(t )  then  has  an  initial  value  of  g/f(0)  =  I  and  a 
value  at  infinity  of  gR(co)  =  G{ oo)  /  G0  (Simulia,  2010). 

In  ABAQUS,  the  viscoelastic  material  is  assumed  to  be  defined  by  a  Prony  series 
expansion  of  the  dimensionless  relaxation  function  as  follows 
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(2.5) 


N 

1=1 

Where  N,  gf ,  and  rf  are  material  constants  with  i=l,  2,  N.  Using  Eq.  (2.5),  the 
expression  for  shear  stress  becomes 

r(<)  =  G0(r-2>,)  (2.6) 

i= 1 

where 

r,=  ^\‘„e"h?r(t-sVs  (2.7) 

The  variables  yt  are  state  variables  that  control  the  stress  relaxation.  The  preceding 

equations  are  used  by  ABAQUS  to  define  the  time-dependent  shear  behavior  in  the 
viscoelastic  model.  To  obtain  the  Prony  series  coefficients  ABAQUS  can  perform  a  curve 
fitting  viscoelastic  material  evaluation  based  on  user  provided  relaxation  data,  as  described 
in  section  2.3.2.  A  similar  method  is  used  for  the  volumetric  response  of  the  system 
(Simulia,  2010). 

II.2.3.2  Calibration  of  Material  Model 

The  numerical  results  in  this  chapter  were  calibrated  by  using  experimental  results 
found  in  literature  for  polymethyl  methacrylate  (PMMA).  All  analyses  were  conducted 
using  ABAQUS  (Simulia,  2010).  In  ABAQUS  a  viscoelastic  material  model  was  used  to 
model  the  PMMA.  Calibration  of  the  material  model  in  ABAQUS  requires  either  the  input 
of  relaxation  data  or  a  set  of  Prony  series  parameters  to  define  the  material  behavior.  For 
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the  viscoelastic  simulations  presented  in  this  Chapter  the  material  models  were  calibrated 
using  test  data  found  in  literature. 

One  set  of  relaxation  data  was  presented  by  Ramesh  Kumar  and  Narasimhan  (2004). 
The  relaxation  response  in  this  paper  was  reported  graphically.  Hence  a  number  of  data 
points  were  picked  from  the  reported  graphs  in  order  to  create  a  data  series  representation 
of  the  relaxation  data.  The  data  was  then  converted  to  a  dimensionless  form  for  the  use  in 
ABAQUS.  Table  II.2-1  shows  the  dimensionless  data  as  extrapolated  from  literature 
(Ramesh  Kumar  and  Narasimhan,  2004).  In  addition  to  the  viscoelastic  data,  an 
instantaneous  Young’s  modulus  of  2600  MPa  and  a  Poisson’s  ratio  of  0.38  are  selected 
based  on  the  infonnation  given  in  the  paper. 


Table  It. 2-1 :  Dimensionless  relaxation  response. 


gR 

1 

0.823 

0.800 

0.785 

0.773 

0.769 

0.767 

0.763 

0.760 

Time 

0.01 

50 

100 

150 

200 

250 

300 

400 

1000 

Viscoelastic  material  response  can  be  defined  through  a  Prony  series  representation. 
ABAQUS  allows  for  evaluation  of  viscoelastic  material  behavior  using  the  dimensionless 
relaxation  test  data.  During  the  material  evaluation,  ABAQUS  perfonns  a  curve  fitting  to 
the  test  data  and  upon  completion  displays  the  Prony  series  parameters  and  material 
response  for  standard  relaxation  and  creep  tests.  The  Prony  series  parameters  calculated  by 
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ABAQUS  for  the  Ramesh  Kumar  and  Narasimhan  relaxation  data  are  shown  in  Table  II.2- 


2. 

Table  II.2-2:  Two-tenn  Prony  series  from  material  analysis. 


i 

g(i)  Prony 

k(i)  Prony 

x(i) 

1 

0.14366 

0 

11.346 

2 

9.83E-02 

0 

115.69 

The  material  model  calibration  was  performed  through  a  simulation  using  a  single  finite 
element  subjected  to  a  uniaxial  displacement  to  create  a  1%  strain  with  the  viscoelastic 
material  model  defined  using  the  aforementioned  Prony  series  parameters.  Using  the  stress 
data  and  the  relation  a(t)=G(t)so  the  relaxation  response  for  the  ABAQUS  model  was 
detennined.  A  comparison  between  the  viscoelastic  model  simulation  and  the  experimental 
data  presented  by  Ramesh  Kumar  and  Narasimhan  is  shown  in  Figure  II.2-1. 
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3000 


Finite  element  viscoelastic  model  simulation 


Figure  II.2-1:  Comparison  of  experimental  relaxation  response  to  viscoelastic  finite 
element  simulation  for  PMMA. 

The  ABAQUS  simulation  exhibited  good  agreement  with  the  experimental  values,  with 
only  slight  differences  noticeable  throughout  the  response.  With  the  simulation  agreeing 
with  the  values  reported  in  literature  the  model  was  determined  to  be  well  calibrated. 

II.2.4  Finite  Element  Modeling 

The  finite  element  nanoindentation  models  in  this  chapter  were  investigated  using 
ABAQUS.  The  models  used  are  2D  axisymmetric  models  using  either  spherical  or  conical 
indenters.  2D  axisymmetric  models  have  the  advantage  of  allowing  quicker  simulations 
than  equivalent  3D  models. 

II.2.4.1  Verification  of  the  Finite  Element  Model 

A  viscoelastic  material  model  of  polymethyl  methacrylate  (PMMA)  was  calibrated 

from  relaxation  data  reported  by  Ramesh  Kumar  and  Narasimhan  (2004)  as  described 
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earlier  in  the  thesis.  Ramesh  Kumar  and  Narasimhan  also  present  results  for  spherical 
indentation.  To  verify  the  2D  axisymmetric  indentation  model  using  ABAQUS,  a  spherical 
indentation  simulation  was  used  to  recreate  the  results  presented  by  Ramesh  Kumar  and 
Narasimhan. 

The  spherical  indentation  model  was  created  by  using  the  axisymmetric  model  from  the 
previous  chapter  and  adapting  it  to  include  a  spherical  indenter.  The  indenter  was  sized  to  a 
radius  of  3.5mm  to  match  that  of  the  experiment  presented  in  Ramesh  Kumar  and 
Narasimhan  (2004)  and  the  specimen  was  sized  to  17.5mm  on  each  side.  The  specimen 
was  modeled  as  a  viscoelastic  isotropic  deformable  solid  with  elastic  parameters  set  to  an 
instantaneous  E=2600MPa,  v=0.38,  and  viscoelastic  behavior  defined  by  the  Prony  series 
parameters  as  shown  in  Table  II. 2-2.  The  indenter  was  modeled  as  an  elastic  isotropic 
defonnable  solid  with  an  elastic  modulus  large  enough  to  prevent  defonnation  of  the 
indenter  during  indentation.  The  indentation  simulation  used  load  control  with  a  triangular 
load  history  with  a  rise  time  of  10s  and  a  maximum  load  of  100N.  Contact  between  the 
indenter  and  specimen  was  modeled  as  frictionless.  The  finite  element  mesh  with  indenter 
is  shown  in  Figure  II.2-2. 
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Figure  II.2-2:  Finite  element  mesh  with  indenter  used  for  spherical  indentation. 


A  viscoelastic  finite  element  analysis  was  performed  on  the  above  described  model. 
From  the  simulation  the  load  and  displacement  values  were  recorded  for  a  direct 
comparison  with  the  results  presented  in  Ramesh  Kumar  and  Narasimhan  (2004).  Two 
different  figures  were  used  for  comparison. 

The  first  figure  used  in  the  comparison  is  a  displacement-time  history  plot.  The 
experimental  results  extrapolated  from  Ramesh  Kumar  and  Narasimhan  (2004)  are 


compared  with  finite  element  results  from  ABAQUS  and  are  shown  in  Figure  II.2-3.  The 
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two  plots  were  in  excellent  agreement  for  a  large  range  of  the  indentation  process.  There 
was  some  slight  deviation  during  the  unloading  portion  of  the  plot  where  the  finite  element 
simulation  experiences  displacements  were  slightly  lower  than  the  experimental  data. 


Figure  II.2-3:  Comparison  of  simulated  and  experimental  spherical  indentation 
displacement-time  responses  for  PMMA. 

A  comparison  of  the  load-displacement  response  was  also  performed.  The  curves  for 
the  experimental  and  finite  element  load-displacement  results  are  show  in  Figure  II.2-4.  In 
this  figure  the  differences  between  the  curves  become  more  noticeable.  There  were  some 
differences  in  both  the  loading  and  unloading  portions  of  the  curves.  The  differences  during 
loading  were  mostly  during  the  initial  stage  of  the  loading  of  load  under  20N,  with  the 
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loads  from  the  experiment  being  slightly  higher.  During  unloading  the  behavior  was  still 
similar  but  there  were  some  more  pronounced  differences.  Particularly  the  recovery  of  the 
material  in  the  experiment  seems  to  lag  behind  the  response  of  the  finite  element 
simulation.  The  latter  difference  is  probably  due  to  a  slight  elastoplastic  behavior  in  the 
material. 


Figure  II.2-4:  Comparison  of  simulated  and  experimental  spherical  indentation  load- 
displacement  responses  for  PMMA. 


In  the  paper  by  Ramesh  Kumar  and  Narasimhan  (2004)  the  experimental  results  were 
similarly  compared  to  numerical  and  analytical  results  for  the  indentation.  The  comparisons 
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performed  in  this  thesis  compare  well  to  those  presented  before  and  similar  differences 
were  observed  between  experimental  and  numerical  results.  The  simulations  presented  here 
then  show  good  verification  of  the  finite  element  model. 
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II.2.5  Results  and  Discussion 


It  has  been  shown  by  Ramesh  Kumar  and  Narasimhan  (2004)  and  Tranchida  et  al. 
(2007)  that  obtaining  viscoelastic  material  parameters  from  nanoindentation  plots  with  the 
conventional  methods  used  for  elastic  and  elastoplastic  materials  has  many  inaccuracies. 
Using  the  Prony  series  parameters  shown  in  Table  II.2-2  as  a  starting  point,  an  investigation 
was  performed  to  see  whether  experimental  nanoindentation  results  on  a  sample  of 
polymethyl  methacrylate  (PMMA)  as  presented  in  a  different  study  could  be  matched 
through  adjustment  of  the  Prony  series  parameters. 

A  study  that  provided  load-displacement  results  from  nanoindentation  on  PMMA  was 
presented  by  Oyen  (2006).  The  nanoindentations  reported  by  Oyen  were  perfonned  by 
using  a  Berkovich  indenter  and  a  triangular  load  history  to  a  maximum  indentation  load  of 
lOmN.  The  experiments  used  for  comparison  had  a  rise  time  of  500s.  To  perfonn  the  finite 
element  simulation  of  the  nanoindentation  a  2D  axisymmetric  model  using  an  equivalent 
70.3°  conical  indenter  was  used.  A  2D  model  was  chosen  due  to  its  lower  computational 
requirements  allowing  for  more  iterations  of  the  simulation  in  a  shorter  time  frame. 

After  running  an  initial  finite  element  analysis  of  the  nanoindentation  using  the  same 
loading  history  as  presented  by  Oyen  (2006)  the  simulation  and  documented  experimental 
results  were  compared.  This  comparison  is  shown  in  Figure  II.2-5.  The  two  response 
curves  were  in  good  agreement  through  the  loading  portion  of  the  indentation,  with  only 
slight  differences  being  noticeable.  However,  during  the  unloading  phase  of  the  indentation 
the  responses  differed  significantly.  The  experiment  shows  a  significant  amount  of  residual 


depth  after  the  completion  of  the  test,  whereas  in  the  simulation  the  residual  depth  was 
fairly  minimal. 


Figure  II.2-5:  Comparison  of  load-displacement  responses  of  PMMA  obtained  from  finite 
element  simulation  against  the  experimental  data  by  Oyen  (2006).  Material  model 
parameters  are  calibrated  from  Ramesh  Kumar  and  Narasimhan  (2004)  data. 


The  simulation  seemed  to  do  a  good  job  of  predicting  the  loading  response  during  the 

indentation;  however  the  unloading  response  did  not  provide  a  good  prediction  of  the 

experimental  behavior.  These  differences  could  have  been  due  to  the  differences  in 

morphology  between  the  PMMA  specimens  used  in  the  two  studies  and/or  due  to  the 
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effects  of  the  sharp  Berkovich  indenter  used  in  obtaining  the  experimental  response.  Using 
the  calibrated  response  as  a  starting  point,  an  investigation  into  the  effect  of  adjusting  the 
Prony  series  parameters  was  performed. 

In  order  to  get  an  understanding  of  how  different  Prony  series  parameters  adjustments 
would  affect  the  load-displacement  response  in  a  nanoindentation  test  a  sensitivity  study 
was  perfonned.  For  the  sensitivity  study  a  single  parameter  was  adjusted  at  a  time  and  the 
change  in  the  shape  of  the  response  curve  was  noted.  Adjustments  were  made  using  a  two- 
term  Prony  series  so  that  adjustments  could  be  compared  to  the  calibrated  response. 


Figure  II.2-6:  Change  in  simulation  response  after  an  increase  in  the  g(2)  term. 
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The  first  adjustment  was  made  to  the  g(2)  term.  The  value  was  adjusted  from  0.0983  to 
0.8.  The  simulation  result  from  this  single  adjustment  is  shown  in  Figure  II.2-6.  The  figure 
shows  that  this  adjustment  produces  a  very  significant  increase  in  both  the  maximum 
indentation  depth  from  the  lOmN  force  as  well  as  the  residual  displacement  after 
unloading. 

An  adjustment  was  then  made  to  the  x(2)  tenn  in  the  Prony  series.  The  term  was 
adjusted  from  the  original  value  of  1 15.69  to  500.  The  difference  in  the  simulation  is  shown 
in  Figure  II.2-7.  This  figure  shows  that  this  adjustment  has  the  effect  of  decreasing  the 
maximum  indentation  depth.  The  two  parameters  have  opposing  effects  in  how  they  affect 
the  simulation  response.  It  was  thought  then  that  adjustments  to  both  parameters  could 
provide  a  material  response  where  the  maximum  indentation  depth  remains  the  same  but 
the  residual  displacement  after  indentation  has  increased.  A  series  of  simulations  were  then 
perfonned  making  adjustments  to  both  parameters  to  test  this. 
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Figure  II.2-7:  Change  in  simulation  response  after  an  increase  in  the  x(2)  term. 


Through  the  series  of  simulations  a  new  finite  element  response  was  obtained  that  was 
closer  to  the  experimental  curve.  Since  making  adjustments  to  increase  the  residual 
displacement  also  had  the  effect  of  increasing  the  maximum  displacement  at  peak  load,  it 
was  difficult  to  maintain  the  same  initial  loading  behavior  which  matched  well  from  the 
initial  calibrated  simulation.  Adjustments  were  ultimately  made  to  both  terms  of  the  Prony 
series  and  the  final  term  values  were;  g(l)=0.1 1,  x(l)=100,  g(2)=0.3,  x(2)=200.  The 
response  using  these  parameter  adjustments  is  shown  in  Figure  II.2-8.  As  the  figure  shows 
the  loading  response  maintained  good  agreement  with  the  original  calibrated  simulation 

and  the  experimental  results,  however,  the  maximum  displacement  had  increased  slightly. 
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The  residual  displacement  had  increased  in  the  new  simulation  to  account  for  50%  of  that 
in  the  experimental  data. 


Figure  II.2-8:  Comparison  of  simulation  before  and  after  Prony  adjustments  to 
experimental  Berkovich  nanoindentation  load-displacement  response  by  Oyen  (2006)  for 
PMMA. 


It  was  shown  that  adjusting  Prony  series  parameters  to  affect  the  shape  of  the  simulated 
load-displacement  curve  to  more  closely  match  an  experimental  curve  was  possible. 
However,  the  adjustments  seemed  to  have  a  limit  to  how  much  they  could  be  adjusted 
while  maintaining  the  behavior  that  was  already  in  good  agreement. 
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After  the  parametric  study,  a  different  relaxation  response  of  PMMA  as  reported  in  a 
paper  by  Sane  and  Knauss  (2001)  was  examined  for  its  response  in  nanoindentation. 
Similar  to  before,  data  points  were  taken  from  a  relaxation  response  figure  and  the  data 
series  was  used  to  define  a  viscoelastic  material.  A  material  analysis  in  ABAQUS  revealed 
that  this  material  was  defined  by  a  six-term  Prony  series.  Figure  II.2-9  shows  the 
indentation  response  of  the  material  in  comparison  to  the  simulation  response  from  the 
Prony  series  parameter  study  and  the  experimental  indentation  curve  from  Oyen  (2006). 
The  loading  response  of  all  considered  simulations  and  the  experimental  data  were  in  good 
agreement.  Interestingly,  the  simulated  response  of  the  PMMA  presented  by  Sane  and 
Knauss  (2001)  was  very  similar  to  the  response  obtained  in  the  Prony  series  parametric 
study,  except  for  some  slight  deviation  during  the  latter  half  of  the  unloading  phase. 
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Figure  II.2-9:  Comparison  of  simulations  from  Prony  series  adjustments  and  that  calibrated 
from  Sane  and  Knauss  (2001)  data  to  experimental  Berkovich  nano  indentation  load- 
displacement  response  for  PMMA. 


The  similarities  in  response  from  the  Prony  series  parametric  study  simulation  and  the 

simulation  calibrated  from  the  data  in  Sane  and  Knauss  (2001)  provided  motivation  to  look 

into  the  relaxation  responses  of  the  finite  element  models  used  in  this  study.  The  relaxation 

response  was  attained  through  finite  element  simulation  as  described  earlier  in  this  thesis.  A 

comparison  of  the  finite  element  analysis  results  for  the  relaxation  responses  of  the 

calibrated  Ramesh  Kumar  and  Narasimhan  (2004)  and  Sane  and  Knauss  (2001)  models  as 

well  as  the  model  developed  through  the  Prony  parameter  study  are  shown  in  Figure  II.2- 

10.  The  comparison  shows  that  the  relaxation  response  for  the  Sane  and  Knauss  data  is  well 

predicted  by  the  model  developed  through  the  Prony  series  parametric  study.  This  showed 
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that  it  is  possible  to  predict  different  morphology  responses  for  the  same  type  of  polymer 
through  Prony  series  parameter  adjustment. 


Figure  II.2-10:  Comparison  of  relaxation  responses  of  PMMA  for  simulations  using  Prony 
adjustments,  calibrated  from  Sane  and  Knauss  (2001)  data,  and  calibrated  from  Ramesh 
Kumar  and  Narasimhan  (2004)  data. 

After  adjusting  Prony  series  parameters  in  an  attempt  to  simulate  an  experimental 
nanoindentation  response  through  numerical  analysis  it  seemed  that  there  was  a  limit  to  the 
extent  of  the  manipulation.  The  adjusted  model  had  a  response  closer  to  the  experimental 
curve  but  there  was  still  a  significant  difference  during  unloading.  However,  the  adjusted 
model  did  match  well  with  the  simulation  of  another  model  calibrated  from  an 
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experimental  relaxation  response  found  in  literature.  It  is  possible  that  for  the  prediction  of 
the  experimental  nanoindentation  curves  a  more  complex  Prony  series  defined  by  more 
terms  would  be  needed.  However,  it  seems  more  likely  that  the  additional  residual 
deformation  during  unloading  is  due  to  the  onset  of  plastic  defonnation  due  to  effects  from 
the  sharp  pyramidal  indenter. 

II.2.6  Conclusions 

A  finite  element  model  was  calibrated  for  polymethyl  methacrylate  (PMMA)  using 
relaxation  data  reported  in  literature.  The  model  displayed  good  agreement  with  the 
reported  relaxation  and  spherical  indentation  responses.  The  nanoindentation  response  of 
the  finite  element  model  using  a  Berkovich  equivalent  conical  indenter  was  then  compared 
to  experimental  results  reported  in  literature  using  a  Berkovich  indenter.  The  comparison 
showed  good  agreement  in  the  loading  response.  However,  the  unloading  response  for  the 
experimental  curve  displayed  a  greater  amount  of  residual  displacement  than  shown  in  the 
numerical  result. 

Prony  series  parameters  for  the  calibrated  model  were  calculated  by  evaluating  the 
material  in  ABAQUS.  A  series  of  Prony  parameters  adjustments  were  perfonned  in  an 
effort  to  more  closely  match  the  numerical  solution  to  the  experimental  nanoindentation 
result.  The  resulting  adjusted  model  had  a  response  more  similar  to  the  experimental  result 
but  still  had  a  smaller  residual  displacement  during  unloading.  It  was  found  however  that 
this  adjusted  model  provided  a  similar  response  to  a  finite  element  model  calibrated  with 
other  PMMA  relaxation  data  reported  in  literature. 
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The  results  show  that  it  is  possible  to  provide  predictions  for  different  samples  of 
PMMA  through  Prony  series  parametric  studies.  However,  there  are  still  issues  in  regard  to 
providing  predictions  of  Berkovich  indenter  nanoindentation.  The  residual  displacements 
experienced  in  the  experimental  indentation  are  thought  to  be  from  plastic  defonnations 
caused  by  the  sharp  pyramidal  geometry  of  the  Berkovich  indenter. 
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II. 3:  Conclusions  and  Future  Challenges 


Elastoplastic  nanoindentations  using  conical  (70.3°)  and  Berkovich  indenters  were 
investigated  using  finite  element  simulations  to  study  the  effect  of  indenter  geometry  on  the 
response  of  the  material.  The  simulations  show  clear  differences  in  the  load-displacement 
responses  of  two  indenters.  Loads  in  the  Berkovich  indentation  are  noticeably  lower  even 
though  the  contact  areas  during  indentation  are  very  similar.  The  result  was  lower 
nonnalized  contact  stress  for  the  Berkovich  indenter.  Examination  of  stress  fields  near  the 
indenter  tip  showed  that  maximum  equivalent  plastic  strain  induced  by  indenter  is  20% 
higher  for  the  Berkovich  indenter.  The  higher  localized  plastic  strains  would  result  in  less 
resistance  to  indentation  by  the  material  and  explain  the  load-displacement  behavior.  The 
results  reveal  that  the  use  of  a  70.3°  conical  indenter  to  approximate  the  response  from  a 
Berkovich  indenter  is  not  valid  for  all  materials. 

The  ability  of  finite  element  analysis  to  predict  experimental  nanoindentation  results 
was  also  examined.  Comparisons  suggest  that  machining-induced  residual  stresses  have 
likely  affected  the  experimental  results.  Residual  stresses  may  be  estimated  by  matching 
simulated  load-displacement  curves  within  the  range  of  experimental  curves. 

Finite  element  analysis  of  viscoelastic  material  was  also  examined.  Using  relaxation 
data  reported  in  literature  a  finite  element  model  was  calibrated  for  polymethyl 
methacrylate  (PMMA).  Simulations  of  a  relaxation  and  spherical  indentation  tests 
displayed  good  agreement  with  the  reported  curves.  Nanoindentation  analysis  using  a 
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Berkovich  equivalent  conical  indenter  showed  good  agreement  in  the  loading  response 
compared  to  reported  experimental  curves  using  a  Berkovich  indenter.  However,  the 
unloading  response  displayed  a  greater  amount  of  residual  displacement  in  the 
experimental  results. 

A  Prony  parameter  adjustment  study  was  perfonned  in  an  effort  to  more  closely  match 
the  numerical  solution  to  the  experimental  nano  indentation  result.  The  adjusted  model  had 
a  response  more  similar  to  the  experimental  result  but  still  had  a  smaller  residual 
displacement  during  unloading.  However,  the  adjusted  model  provided  a  similar  response 
to  a  finite  element  model  calibrated  with  different  PMMA  relaxation  data  reported  in 
literature. 

It  was  shown  to  be  possible  to  provide  predictions  for  different  samples  of  PMMA 
through  Prony  series  parametric  studies.  Conversely,  there  are  challenges  in  providing 
predictions  of  Berkovich  indenter  nanoindentation.  The  residual  displacements  in  the 
experimental  indentation  are  considered  to  be  from  plastic  deformations  caused  by  the 
sharp  pyramidal  geometry  of  the  Berkovich  indenter. 

Insights  have  been  gained  into  the  behavior  of  elastoplastic  and  viscoelastic  materials 
during  nanoindentation.  Nevertheless,  there  are  many  questions  that  remain  in  the  areas  of 
numerical  and  experimental  nanoindentation.  This  thesis  has  explored  the  comparison  of 
conical  and  Berkovich  indenters  through  simulation.  The  examination  of  additional 
materials  following  various  hardening  laws  is  still  needed.  In  addition,  investigations  into 
the  comparison  of  the  two  “equivalent”  tip  geometries  through  experimental 
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nanoindentations  in  different  materials  would  provide  valuable  infonnation  into  the 
equivalence  assumption. 

In  the  viscoelastic  comparisons,  plastic  deformation  effects  were  observed  in  the 
experimental  indentations.  Further  work  examining  the  effect  of  plastic  defonnation  in 
polymer  indentations  is  needed.  Methods  for  developing  a  finite  element  viscoplasticity 
model  based  on  experimental  nanoindentation  curves  would  be  particularly  valuable. 

The  methods  used  for  the  examination  of  material  behavior  in  nanoindentation  have 
been  extensively  developed  since  the  inception  of  the  technique.  The  methods  work  well 
for  many  types  of  materials  and  for  different  indenter  types.  Modern  nanoindenters  and 
finite  element  software  are  making  examination  of  these  problems  more  accessible. 
Flowever,  there  are  still  many  challenges  that  remain  for  the  examination  of  more  complex 
material  types  and  how  they  will  respond  to  different  indenter  geometries. 
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